C 20 AUG 07 - LVS - PBC-EFP CLEANED
C 24 MAR 07 - LVS - PAD FRGINF,EFPPAR; PRINT FUNCTION PREFP ADDED
C 17 JAN 06 - HL  - MD SWITCHING FUNCTIONS, CHANGES TO IND.DIP.,EFP/PCM
C 14 NOV 05 - HL  - POLPCM ADDED TO ALLOW FOR EFP+PCM COMPUTATIONS
C 19 SEP 05 - HMN,IA - CHANGES FOR MOLECULAR DYNAMICS
C 15 JUL 05 - MWS - EFPPOL,EFPPOL2: RESTRICT VCLR RANGES
C 13 JUL 05 - MWS - CHANGE DYN.STORE FOR CHGTRN,EFPPOL
C  5 JUL 05 - MWS - ALLOCATE DYNAMIC MEMORY FOR DIPIT
C 27 JUN 05 - HL  - PERIODIC BOUNDARY COND, PARALLEL EFP2 DISPERSION
C  1 JUN 05 - MWS - EFSP: SUM ALL ES TERMS IN PARALLEL RUN
C 30 APR 05 - DGF - SYNCHRONISE OPTGRD COMMON BLOCK
C 14 MAR 05 - HL  - EFP CHARGE TRANSFER ENERGY, PAULI CHANGES
C  7 MAR 05 - IA  - EFP2 DISPERSION AND DYNAMIC POLARIZABILITY
C 22 FEB 05 - MWS - EFSP: CHANGE ENERGY PRINTING (AND DEBUG OUTPUTS)
C 13 FEB 05 - HMN - ADD DUMMY EFP MD HOOKS INTO EFSP,
C                   ADD EFP PARALLELIZATION INTO EFSP AND EFPPOL
C 13 FEB 05 - JS  - ADD EFP MP2 CALLS TO EFSP AND EFPPOL, ADD EFPPOL2
C 23 JUL 04 - HL  - EFSP: DO FORMAL MEMORY ALLOCATION FOR EFPPOL CALL
C  9 DEC 03 - MWS - EFGRAD: KILL NUMERICAL HESSIAN REQUEST
C  3 JUL 03 - JMM - SUPPRESS PRINTING FOR MONTE CARLO JOBS
C  7 AUG 02 - HT  - USE DYNAMIC MEMORY FOR EFP PAULI REPULSION
C 17 APR 02 - MWS - SYNCH UP FRGINF COMMON
C 13 JUN 01 - MAF - EFSP: COMPUTE CHARGE PENETRATION ENERGY
C 20 FEB 01 - PND - EFSP,EFPPOL: DIP.CONV. AND PRINTING CHANGES
C 25 MAR 00 - MWS - CHANGE ARGUMENTS TO VNNDER
C 20 NOV 97 - MWS - EFSP: USE JAN'S PAULI EXCHANGE ROUTINE
C  2 SEP 97 - MWS - EFSP: CHANGE WORDING OF ENERGY OUTPUT
C 16 MAY 97 - GNM - EFFECTIVE FRAGMENT DRIVER MODULE CREATED
C
C*MODULE EFDRVR  *DECK EFSP
      SUBROUTINE EFSP
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      LOGICAL IFGRAD
      LOGICAL IFEWLD,TNFOIL,EWLDPL,DEBUGOPT
C
      PARAMETER (MXAO=8192, MXATM=2000, MXDFG=5, MXFRG=1050,
     *           MXPT=2000, MXFGPT=12000, MXIFRQ=12,
     *           MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO = 0.0D+00, HALF=0.50D+00, FACT=1.593601268D-03)
C
      DOUBLE PRECISION MD
C
      COMMON /EFPCMO/ NOCC(MXDFG),NVIR(MXDFG),NAO(MXDFG),NTAO,
     *                CANFOK(MXAO,MXDFG)
      COMMON /EFPIO / MULMAT(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 /FRGNRG/ CCNRG,CDNRG,CQNRG,CONRG,DDNRG,DQNRG,QQNRG,
     *                ESNRG,REPNRG,POLNRG,POLNRG2,DISNRG,DISNRG8,
     *                EDISD,CHTNRG
      COMMON /FUNCT / E,EG(3*MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OPTGRD/ Q(3*MXATM),ES,QQ(20),
     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGH
      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
      COMMON /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
     *                NPRTGO,IDPUNC,IGOFLG
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,KEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IRPPCM,IEF,IP_F,NFMOPCM
      COMMON /PCMPNT/ LXYZRE,LSSFE,LLIST,LALPCM,LRINPCM,LINAPCM,LINFPCM,
     *                LINIPCM,LMEPCM,LAXYZCT,LXYZCT2,LNVERT,LQSN,LQSE,
     *                LQFS,LQIND,LISPHE,LVAD,LQOR,LVECMUL,LAIPRJ,LFIPRJ
     *                ,LPEL
      COMMON /PCMPRT/ GCAVP,GCAVS,GDISP,GREP,EHFGAS
      COMMON /PCMEFP/ PB_MUL,PC_MUL,PB_POL,PC_POL,P_FF,P_NF,P_FN,P_I,
     *                P_IBIS,P_J,P_NUCC,P_NUCCBIS,ENPCM
      COMMON /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 /EWALD / IFEWLD,TNFOIL,BETA,ELCNST,EPSLN,CUTOFF,CUTLST,
     *                KMAX,NKVEC,LEVEL,EWLDPL
      COMMON /SELFE / ESELFCC, ESELFDD, EEXCLCC, EEXCLCD, EEXCLDD,
     *                EEXCLCQ 
      COMMON /EFPMP2/ MP2TYP
C
      DATA HSSIAN/8HHESSIAN /, SADPT/8HSADPOINT/
      DATA AIRC/8HIRC     /
      DATA MD /8HMD      /
C
C     -------------------------------------------------------------
C     --- THIS ROUTINE CALCULATES FRAGMENT-FRAGMENT INTERACTION ---
C     --- ENERGIES IN THE ABSENCE OF AN AB INITIO MOLECULE      ---
C     -------------------------------------------------------------
C
C     CCNRG  = CHARGE-CHARGE INTERACTION ENERGY
C     CDNRG  = CHARGE-DIPOLE INTERACTION ENERGY
C     CQNRG  = CHARGE-QUADRUPOLE INTERACTION ENERGY
C     CONRG  = CHARGE-OCTUPOLE INTERACTION ENERGY
C     DDNRG  = DIPOLE-DIPOLE INTERACTION ENERGY
C     DQNRG  = DIPOLE-QUADRUPOLE INTERACTION ENERGY
C     QQNRG  = QUADRUPOLE-QUADRUPOLE INTERACTION ENERGY
C     ESNRG  = TOTAL ELECTROSTATIC INTERACTION ENERGY
C     REPNRG = EXCHANGE-REPULSION ENERGY
C     POLNRG = POLARIZATION ENERGY
C     POLNRG2= MP2/EFP1 2ND ORDER POLARIZATION ENERGY CORRECTION
C     DISNRG = MP2/EFP1 C6 DISPERSION ENERGY
C     DISNRG8= MP2/EFP1 C8 DISPERSION ENERGY
C     EDISD  = EFP2 DAMPED DISPERSION ENERGY
C     CHTNRG = EFP2 CHARGE TRANSFER ENERGY
C     ENPCM  = PCM ELECTROSTATIC SOLVATION ENERGY
C     GCAVP  = PCM CAVITATION ENERGY
C     GDISP  = PCM DISPERSION ENERGY
C     GREP   = PCM REPULSION ENERGY
C
      E      = ZERO
      CCNRG  = ZERO
      CDNRG  = ZERO
      CQNRG  = ZERO
      CONRG  = ZERO
      DDNRG  = ZERO
      DQNRG  = ZERO
      QQNRG  = ZERO
      ESNRG  = ZERO
      REPNRG = ZERO
      POLNRG = ZERO
      POLNRG2= ZERO
      DISNRG = ZERO
      DISNRG8= ZERO
      EDISD  = ZERO
      CHTNRG = ZERO
      ENPCM  = ZERO
      GCAVP  = ZERO
      GDISP  = ZERO
      GREP   = ZERO
C
C     -- EFP CENTER OF MASS --
C
      CALL EFCM
C
C     -- ELECTROSTATIC ENERGIES --
C     ELECTROSTATIC ENERGIES (AND GRADIENTS) ARE DONE IN FFELEC NOW,
C     REPLACING THE SMALL ARMY OF ROUTINES CHGCHG, ETC. ETC.
C
C      WRITE(IW,*)'DESFRG-0 ',(DESFRG(I,1), I=1,6)
      MAXDER = 0
      IFGRAD = .FALSE.
      CALL DERCHK(MAXDER)
      IF (MAXDER.GT.0) IFGRAD = .TRUE.
      CALL FFELEC(ESNRG,CCNRG,CDNRG,CQNRG,CONRG,DDNRG,DQNRG,QQNRG,
     *                  IFGRAD,.FALSE.)
C
      IF (IFEWLD) THEN
C  ---- EWALD SUMS: DIRECT TERM IS CALCULATED WITHIN FFELEC
C  ---- HERE WE CALCULATE OTHER TERMS:
C  ---- SELF-ENERGY & EXCLUSION ENERGY, RECIPROCAL, AND TIN-FOIL PART
C
C      WRITE(IW,*)'DESFRG-1 ',(DESFRG(I,1), I=1,6)
C      WRITE(6,*)'ESNRG-1 = ',ESNRG
C
C         CALL SELFCC(ESELFCC,ESELFDD,EEXCLCC,EEXCLCD,EEXCLDD,IFGRAD)
         ESNRG = ESNRG + ESELFCC + ESELFDD + EEXCLCC + EEXCLCD 
     *         + EEXCLDD + EEXCLCQ 
 
C      WRITE(IW,*)'DESFRG-2 ',(DESFRG(I,1), I=1,6)
C      WRITE(6,*)'ESNRG-2 = ',ESNRG
C
C         CALL KVECTORS
C
C         CALL VALFM(LOADFM)
C         WRITE(6,*)'MEMORY AFTER KVECTORS=',LOADFM
C         LCOSCH = LOADFM + 1
C         LSINCH = LCOSCH + MAXVEC
C         LTCHCH = LSINCH + MAXVEC
C         LAST  = LTCHCH + MAXVEC
C         NEED = LAST - LOADFM - 1
C         WRITE(6,*)'MEMORY FOR PECIPR NEEDED=',NEED
C         CALL GETFM(NEED)
C         WRITE(6,*)'MEMORY FOR PECIPR ALLOCATED'
         CALL RECIPR(ERCPRCC,ERCPRCD,ERCPRDD,ERCPRCQ,ERCPRDQ,
     *                    ERCPRQQ,ERCPRCO,IFGRAD)
C         CALL RETFM(NEED)
C         WRITE(6,988) ERCPRCC,ERCPRCD,ERCPRDD,ERCPRCQ
C  988  FORMAT(/1X,'ERCPRCC = ',F18.14/
C     *        1X,'ERCPRCD = ',F18.14/
C     *        1X,'ERCPRDD = ',F18.14/
C     *        1X,'ERCPRCQ = ',F18.14)
         ESNRG = ESNRG + 
     *           ERCPRCC + ERCPRCD + ERCPRDD  + ERCPRCQ 
C
C      IF (MASWRK) THEN
C      WRITE(IW,*)'DESFRG-3 ',(DESFRG(I,1), I=1,6)
C      WRITE(6,*)'ESNRG-3 = ',ESNRG
C      END IF
        CALL TINFOIL(ETFCC,ETFCD,ETFDD, IFGRAD)
C        WRITE(6,989) ETFCC,ETFCD,ETFDD
C  989  FORMAT(/1X,'ETFCC = ',F18.14/
C     *        1X,'ETFCD = ',F18.14/
C     *        1X,'ETFDD = ',F18.14)
         ESNRG = ESNRG + ETFCC + ETFCD + ETFDD      
C       IF (MASWRK) THEN
C         WRITE(IW,*)'DESFRG-FINAL ',(DESFRG(I,1), I=1,6)
C         WRITE(6,*)'ESNRG-FINAL = ',ESNRG
C       END IF
C         WRITE(6,*) 'TOTAL CHARGE-CHARGE: ',CCNRG + ESELFCC + 
C     *               EEXCLCC + ERCPRCC + ETINFOIL
C  ---- THIS IS DEBUGGING OPTION WHICH CALCULATES ELECTROSTATIC SUM 
C  ---- DIRECTLY. FOR UNCLEAR REASON, WORKS PROPERLY ONLY FOR ATOMS NOT 
C  ---- FOR MOLECULES
        DEBUGOPT=.FALSE.
        IF (DEBUGOPT) CALL DIRECTSUM
      END IF
C
C      CALL CHGCHG(CCNRG)
C      CALL CHGDIP(CDNRG)
C      CALL CHGQUA(CQNRG)
C      CALL CHGOCT(CONRG)
C      CALL DPLDPL(DDNRG)
C      CALL DPQUAD(DQNRG)
C      CALL QUDQUD(QQNRG)
C
C      IF(GOPARR) THEN
C        ESBUFF(1) = CCNRG
C        ESBUFF(2) = CDNRG
C        ESBUFF(3) = CQNRG
C        ESBUFF(4) = CONRG
C        ESBUFF(5) = DDNRG
C        ESBUFF(6) = DQNRG
C        ESBUFF(7) = QQNRG
C        CALL DDI_GSUMF(666,ESBUFF,7)
C        CCNRG = ESBUFF(1)
C        CDNRG = ESBUFF(2)
C        CQNRG = ESBUFF(3)
C        CONRG = ESBUFF(4)
C        DDNRG = ESBUFF(5)
C        DQNRG = ESBUFF(6)
C        QQNRG = ESBUFF(7)
C      ENDIF
C
C     -- EXCHANGE-REPULSION ENERGY ---
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
C
         CALL VALFM(LOADFM)
         LAREL = LOADFM + 1
         LCREL = LAREL + LEN
         LEFLD = LCREL + LEN
         LEFAD = LEFLD + 3*NPTTPT
         LAST  = LEFAD + 3*NPTTPT
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
         CALL DAREAD(IDAF,IODA,XX(LAREL),2*LEN,75,0)
         CALL EREPUL(REPNRG,XX(LAREL),XX(LCREL),MXRPTS,NDFRG2)
         CALL RETFM(NEED)
      ELSE IF (NTMO.GT.0) THEN
         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
C
         IF(RUNTYP.NE.MD) THEN
           CALL VALFM(LOADFM)
           LPROVEC = LOADFM  + 1
           LFOCKMA = LPROVEC + MXBF*NTMO
           LSMAT   = LFOCKMA + MXMO2*NFRG
           LTMAT   = LSMAT   + MXBF*MXBF
           LWRK    = LTMAT   + MXBF*MXBF
           LSIJ    = LWRK    + MXBF
           LTIJ    = LSIJ    + MXMO*MXMO
           LFASQ   = LTIJ    + MXMO*MXMO
           LFBSQ   = LFASQ   + MXMO*MXMO
           LAST    = LFBSQ   + MXMO*MXMO
           NEED=LAST-LOADFM-1
           CALL GETFM(NEED)
           CALL DAREAD(IDAF,IODA,XX(LPROVEC),LENPV,267,0)
           CALL DAREAD(IDAF,IODA,XX(LFOCKMA),LENFM,268,0)
           CALL EXREP(REPNRG,XX(LPROVEC),XX(LFOCKMA),XX(LSMAT),
     *                XX(LTMAT),XX(LWRK),XX(LSIJ),XX(LTIJ),
     *                XX(LFASQ),XX(LFBSQ),MXBF,MXMO)
           CALL RETFM(NEED)
         ELSE
C          - IF MD, RUN GRADIENTS HERE
           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)
           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(REPNRG,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)
         END IF
C
      END IF
C
C     -- MP2/EFP1 DISPERSION ENERGY --
C
      IF(IEFP2.EQ.1 .OR. MP2TYP.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
         DISNRG = ZERO
         DISNRG8= ZERO
C
         CALL VALFM(LOADFM)
         LAREL = LOADFM + 1
         LDREL = LAREL + LEN
         LCREL = LDREL + LEN
         LAST  = LCREL + LEN
         LBREL = LAST  + LEN
         LEREL = LBREL + LEN
         LAST1 = LEREL + LEN
         NEED  = LAST1 - LOADFM - 1
         CALL GETFM(NEED)
         CALL DAREAD(IDAF,IODA,XX(LAREL),6*LEN,81,0)
         CN=6.0D+00
         CALL EDISPR(DISNRG,XX(LAREL),XX(LDREL),XX(LCREL),
     *               MXDPTS,NDFRG2,CN)
         CN=8.0D+00
         CALL EDISPR(DISNRG8,XX(LAST),XX(LBREL),XX(LEREL),
     *               MXDPTS,NDFRG2,CN)
         CALL RETFM(NEED)
      END IF
C
C     -- POLARIZATION ENERGY ---
C     -- PCM SOLVATION ENERGY --
C
      IF(IEFP2.EQ.1) THEN
        CALL VALFM(LOADFM)
        LEFLD = LOADFM + 1
        LEFAD = LEFLD + 3*NPTTPT
        LAST  = LEFAD + 3*NPTTPT
        NEED=LAST-LOADFM-1
        CALL GETFM(NEED)
        IEFPOL=1
        CALL EFPPOL(XX(LEFLD),XX(LEFAD),POLNRG)
        IEFPOL=0
        IEFPOL2=1
        CALL EFPPOL2(XX(LEFLD),XX(LEFAD),POLNRG2)
        IEFPOL2=0
        CALL RETFM(NEED)
      ELSE
        IF(IPCM.EQ.1 .OR. IPCM.EQ.817) THEN
          IPCM = 817
          NPRTBK=NPRINT
          NPRINT=817
          IF(ICAV.EQ.1) CALL CAVITM
          IF(IDISP.EQ.1) CALL DISRPM
          NPRINT=NPRTBK
C
          CALL VALFM(LOADFM)
          LEFLD = LOADFM + 1
          LEFAD = LEFLD  + 3*NPTTPT
          LAST  = LEFAD  + 3*NPTTPT
          NEED=LAST - LOADFM -1
          CALL GETFM(NEED)
          CALL POLPCM(XX(LEFLD),XX(LEFAD),XX(LQSE),
     *                XX(LQSN),XX(LAXYZCT),XX(LAXYZCT+MXTS),
     *                XX(LAXYZCT+MXTS*2),XX(LAXYZCT+MXTS*3),
     *                XX(LXYZRE),XX(LXYZRE+MXSP),XX(LXYZRE+MXSP*2),
     *                XX(LXYZRE+MXSP*3),XX(LISPHE),XX(LLIST),POLNRG)
          CALL RETFM(NEED)
        ELSE
          CALL VALFM(LOADFM)
          LEFLD = LOADFM + 1
          LEFAD = LEFLD + 3*NPTTPT
          LAST  = LEFAD + 3*NPTTPT
          NEED=LAST-LOADFM-1
          CALL GETFM(NEED)
          CALL EFPPOL(XX(LEFLD),XX(LEFAD),POLNRG)
          CALL RETFM(NEED)
        END IF
      END IF
C
C     -- EFP2 CHARGE TRANSFER ENERGY --
C
C EFP2 CHARGE TRANSFER NOT YET DERIVED FOR OPEN SHELL (ROHF) METHOD
C
C CHECK TO SEE IF ANY MULTIPLICITIES ARE GREATER THAN 1
C IF SO, SKIP CT CALCULATIONS ALTOGETHER
C
C THIS SKIPS ANY CT CALCULATION BETWEEN
C CLOSED-SHELL MOLECULES AS WELL
C
      DO I=1,NFRG
       IF (MULMAT(I) .GT. 1)  GOTO 100
      END DO
C
      IF(NTAO.GT.0)THEN
         MXAOEFP=0
         MXBFEFP=0
         DO I = 1, NFRG
            MXAOEFP=MAX(MXAOEFP,NAO(ISET(I)))
            MXBFEFP=MAX(MXBFEFP,NPBF(I))
         END DO
         LENCV=MXBFEFP*NTAO
C
         IF(RUNTYP.NE.MD) THEN
           CALL VALFM(LOADFM)
           LCANVEC = LOADFM  + 1
           LSMAT   = LCANVEC + MXBFEFP*NTAO
           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
           LECTFRG = LVBBEFA + MXAOEFP*MXAOEFP
           LIDOIT  = LECTFRG + NFRG*NFRG
           LQUAD   = LIDOIT  + NMTTPT
           LAST    = LQUAD   + 6*NMTTPT
           NEED=LAST-LOADFM-1
           CALL GETFM(NEED)
           CALL DAREAD(IDAF,IODA,XX(LCANVEC),LENCV,269,0)
           CALL CHGTRN(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(LECTFRG),
     *                 XX(LIDOIT),XX(LQUAD),MXBFEFP,MXAOEFP)
           CALL RETFM(NEED)
         ELSE
C          - IF MD, RUN GRADIENTS HERE
           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)
           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)
         END IF
      END IF
C
C     -- EFP2 DYNAMIC POLARIZATION ENERGY/GRADIENTS ---
C
C DDK
C EFP2 DISPERSION NOT YET DERIVED FOR  OPEN SHELL (ROHF) METHOD
C
C CHECK TO SEE IF ANY MULTIPLICITIES ARE GREATER THAN 1
C IF SO, SKIP DISPERSION CALCULATIONS ALTOGETHER
C
C THIS SKIPS ANY DISPERSION CALCULATION BETWEEN
C CLOSED-SHELL MOLECULES AS WELL
C
      DO I=1,NFRG
       IF (MULMAT(I) .GT. 1)  GOTO 100
      END DO
      IF (NDPTTPT.GT.0) THEN
         MXPPTS = 0
         DO IFRG = 1,NFRG
            MXPPTS = MAX(NPPTS(IFRG),MXPPTS)
         ENDDO
         CALL VALFM(LOADFM)
         LDYNTR = LOADFM + 1
         LAST   = LDYNTR + NFRG*MXIFRQ*MXPPTS
         NEED   = LAST -LOADFM -1
         CALL GETFM(NEED)
         CALL EFPDYN(EDISD,XX(LDYNTR),MXPPTS)
         CALL RETFM(NEED)
      END IF
C
 100  CONTINUE
C
C     ******************
C     -- TOTAL ENERGY --
C     ******************
C         ESNRG = CCNRG + CDNRG + CQNRG + CONRG + DDNRG + DQNRG + QQNRG
C
      E = ESNRG+REPNRG+POLNRG+EDISD+CHTNRG+POLNRG2+DISNRG+DISNRG8
     *  + HALF*ENPCM + FACT*(GCAVP+GDISP+GREP)
C
C
      NEVALS = NEVALS + 1
C
C     -- PRINT OUT ENERGIES --
C
      IF(RUNTYP.EQ.MD) GOTO 800
C
C         IF YOU ASSIGN IGOFLG=2 HERE, YOU WILL GET FULL ENERGY OUTPUT
C         FOR ALL GEOMETRIES COMPUTED, GOOD FOR DEBUGGING TERM BY TERM.
C
C      IGOFLG=2
C
      IF(MASWRK  .AND.  IGOFLG.EQ.2) WRITE(IW,9010)
C
C         NOTE THAT THE TOTAL ENERGY IS INSIDE THE IGOFLG=2 PORTION
C         OF THE ELSE CLAUSE, TO SUPPRESS MONTE CARLO ENERGIES,
C         BUT THE MORE BASIC RUNS ALWAYS PRODUCE TOTAL ENERGY OUTPUT.
C
C     IF((RUNTYP.EQ.OPTIMIZ .AND. NSERCH.NE.0) .OR.
      IF(
     *   (RUNTYP.EQ.SADPT   .AND. NSERCH.NE.0) .OR.
     *   (RUNTYP.EQ.AIRC)                      .OR.
     *   (RUNTYP.EQ.HSSIAN)) THEN
         IF(MASWRK  .AND.  IGOFLG.EQ.2) THEN
            WRITE(IW,9090) ESNRG
            WRITE(IW,9100) REPNRG
            WRITE(IW,9110) POLNRG
            IF(IEFDP.EQ.1) WRITE(IW,9111) EDISD
            IF(NTAO.GT.0) WRITE(IW,9112) CHTNRG
            IF(IEFP2.EQ.1) WRITE(IW,9120) POLNRG2
            IF(IEFP2.EQ.1 .OR. MP2TYP.EQ.1) THEN
              WRITE(IW,9130) DISNRG
              WRITE(IW,9140) DISNRG8
            END IF
            WRITE(IW,9210)
         END IF
         IF(MASWRK) THEN
           IF(IPCM.EQ.1 .OR. IPCM.EQ.817) THEN
              WRITE(IW,9300) E
           ELSE
              WRITE(IW,9200) E
           END IF
         END IF
         IGOFLG=1
      ELSE
         IF(MASWRK.AND.  IGOFLG.EQ.2) THEN
            WRITE(IW,9020) CCNRG
            WRITE(IW,9030) CDNRG
            WRITE(IW,9040) CQNRG
            WRITE(IW,9050) CONRG
            WRITE(IW,9060) DDNRG
            WRITE(IW,9070) DQNRG
            WRITE(IW,9080) QQNRG
            WRITE(IW,9220)
            WRITE(IW,9090) ESNRG
            WRITE(IW,9100) REPNRG
            WRITE(IW,9110) POLNRG
            IF(IEFDP.EQ.1) WRITE(IW,9111) EDISD
            IF(NTAO.GT.0) WRITE(IW,9112) CHTNRG
            IF(IEFP2.EQ.1) WRITE(IW,9120) POLNRG2
            IF(IEFP2.EQ.1 .OR. MP2TYP.EQ.1) THEN
              WRITE(IW,9130) DISNRG
              WRITE(IW,9140) DISNRG8
            END IF
            WRITE(IW,9210)
            IF(IPCM.EQ.1 .OR. IPCM.EQ.817) THEN
               WRITE(IW,9300) E
            ELSE
               WRITE(IW,9200) E
            END IF
         END IF
C        IGOFLG=1
      END IF
C
 800  CONTINUE
C
      RETURN
C
 9010 FORMAT(/17X,38('-')/
     *       17X,'FRAGMENT-FRAGMENT INTERACTION ENERGIES'/
     *       17X,38('-')/)
 9020 FORMAT(17X,'CHARGE-CHARGE         = ',F14.10)
 9030 FORMAT(17X,'CHARGE-DIPOLE         = ',F14.10)
 9040 FORMAT(17X,'CHARGE-QUADRUPOLE     = ',F14.10)
 9050 FORMAT(17X,'CHARGE-OCTUPOLE       = ',F14.10)
 9060 FORMAT(17X,'DIPOLE-DIPOLE         = ',F14.10)
 9070 FORMAT(17X,'DIPOLE-QUADRUPOLE     = ',F14.10)
 9080 FORMAT(17X,'QUADRUPOLE-QUADRUPOLE = ',F14.10)
 9090 FORMAT(17X,'ELECTROSTATIC ENERGY  = ',F14.10)
 9100 FORMAT(17X,'REPULSION ENERGY      = ',F14.10)
 9110 FORMAT(17X,'POLARIZATION ENERGY   = ',F14.10)
 9111 FORMAT(17X,'DISPERSION ENERGY     = ',F14.10)
 9112 FORMAT(17X,'CHARGE TRANSFER ENRGY = ',F14.10)
 9120 FORMAT(17X,'MP2 POLARIZATION ENGY = ',F14.10)
 9130 FORMAT(17X,'C6 DISPERSION ENERGY  = ',F14.10)
 9140 FORMAT(17X,'C8 DISPERSION ENERGY  = ',F14.10)
 9200 FORMAT(17X,'FINAL EFP ENERGY      = ',F14.10)
 9210 FORMAT(17X,38('-'))
 9220 FORMAT(41X,14('-'))
 9300 FORMAT(17X,'FINAL EFPPCM ENERGY   = ',F14.10)
C
      END
C*MODULE EFDRVR  *DECK EFPPOL
      SUBROUTINE EFPPOL(EFLD,EFADD,POLNRG)
C
C --------------------------------------------------------------
C --- THIS ROUTINE CALCULATES FRAGMENT-FRAGMENT POLARIZATION ---
C ---   ENERGIES IN THE ABSENCE OF AN AB INITO MOLECULE      ---
C --------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXPT=2000, MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO=0.0D+00, DEBYE=2.541766D+00)
C
      CHARACTER*8 POLNAM,FRGNME,DNAME,DPSNAM,DPOLNAM
C
      LOGICAL MINMEM
      LOGICAL GOPARR,DSKWRK,MASWRK
      LOGICAL IFEWLD,TNFOIL,EWLDPL,DEBUGOPT
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 /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 /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INDDIP/ TMUX,TMUY,TMUZ,TMUXD,TMUYD,TMUZD,MINMEM
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
     *                NPRTGO,IDPUNC,IGOFLG
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
C
      COMMON /EWALD / IFEWLD,TNFOIL,BETA,ELCNST,EPSLN,CUTOFF,CUTLST,
     *                KMAX,NKVEC,LEVEL,EWLDPL
      COMMON /EWFLD / RCPFLD(3,MXFGPT),RCPIIFLD1(3,MXFGPT),
     *                RCPIIFLD2(3,MXFGPT), TNFLD(3),
     *                TNIIFLD1(3),TNIIFLD2(3),
     *                SELFFLD1(3,MXFGPT),SELFFLD2(3,MXFGPT),
     *                EXCLFLD(3,MXFGPT),
     *                EXCLIIFLD1(3,MXFGPT),EXCLIIFLD2(3,MXFGPT)
C
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
C
      COMMON /POL2  / CHGHF(2,MXFGPT),DIPHF(3,MXFGPT),QADHF(6,MXFGPT),
     *                AMULTHF(MXFGPT),BMULTHF(MXFGPT),
     *                AMULT2HF(MXFGPT),BMULT2HF(MXFGPT),
     *                EFPOL2(9,MXFRG),DIND2(3,MXFGPT),DINDD2(3,MXFGPT),
     *                DCORD(3,MXPT),CDIS(3,MXFGPT),
     *                CLPD(4*MXFGPT),ZLPD(4*MXFGPT),
     *                NLPD(4*MXFGPT),KFD(MXFGPT),KLD(MXFGPT),
     *                DNAME(MXPT),DPSNAM(MXFGPT)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /EFPAR / LOCNFG,NLOW,NUP
C
      DIMENSION CHGTMP(2,MXFGPT),DIPTMP(3,MXFGPT),QADTMP(6,MXFGPT),
     *          EFATMP(MXFGPT),EFBTMP(MXFGPT)
C
      DIMENSION EFLD(3,NPTTPT),EFADD(3,NPTTPT)
C
      DATA ENERGY  /8HENERGY  /,GRADIENT/8HGRADIENT/
C
C     ATTENTION: THIS ROUTINE WILL CALCULATE EPOL(HF) IN MP2/EFP1.
C     HOWEVER, WE SHOULD BE AWARE THAT EFPOL IS HF VALUE BUT EFQAD
C     (IN FRGFLD),EFATRM/EFBTRM(IN FRGFLD),EFCHG(IN FRGFLD),AND EFDIP
C     (IN FRGFLD) ARE MP2 DATA. SO THE CORRESPONDING HF DATA SHOULD BE
C     WRITTEN TO COMMON BLOCK AT THE BEGINNING OF THIS ROUTINE AND
C     AT THE END THE ORIGINAL MP2 DATA SHOULD BE RESTORED.
C
C     BE VERY CAREFUL ABOUT THE DATA UNDER COMMON FRFD, INDDIP,
C     AND EFPPAR (DIND AND DINDD), SINCE THEY ARE WRITTEN WHEN EPOL
C     IS CALCULATED!!!
C
C
      IF((IEFP2.EQ.1).AND.(IEFPOL.EQ.1)) THEN
        DO I=1,MXFGPT
C         SAVE MP2 DATA TO TMP
          EFATMP(I)=EFATRM(I)
          EFBTMP(I)=EFBTRM(I)
          CHGTMP(1,I)=EFCHG(1,I)
          CHGTMP(2,I)=EFCHG(2,I)
          DIPTMP(1,I)=EFDIP(1,I)
          DIPTMP(2,I)=EFDIP(2,I)
          DIPTMP(3,I)=EFDIP(3,I)
          QADTMP(1,I)=EFQAD(1,I)
          QADTMP(2,I)=EFQAD(2,I)
          QADTMP(3,I)=EFQAD(3,I)
          QADTMP(4,I)=EFQAD(4,I)
          QADTMP(5,I)=EFQAD(5,I)
          QADTMP(6,I)=EFQAD(6,I)
C         WRITE HF DATA TO COMMON BLOCK
          EFATRM(I)=AMULTHF(I)
          EFBTRM(I)=BMULTHF(I)
          EFCHG(1,I)=CHGHF(1,I)
          EFCHG(2,I)=CHGHF(2,I)
          EFDIP(1,I)=DIPHF(1,I)
          EFDIP(2,I)=DIPHF(2,I)
          EFDIP(3,I)=DIPHF(3,I)
          EFQAD(1,I)=QADHF(1,I)
          EFQAD(2,I)=QADHF(2,I)
          EFQAD(3,I)=QADHF(3,I)
          EFQAD(4,I)=QADHF(4,I)
          EFQAD(5,I)=QADHF(5,I)
          EFQAD(6,I)=QADHF(6,I)
        ENDDO
      END IF
C
      CALL FRGFLD
      IF (IFEWLD .AND. EWLDPL) THEN 
C   ----  DEBUG: CALCULATES STATIC FIELD EXPLICITELY      
         DEBUGOPT=.FALSE.
         IF (DEBUGOPT) CALL DIRECTSUM_POL
C   ---- TINFOIL CONTRIBUTION TO THE FIELD IS CALCULATED IN TNFOIL, 
C   ---- TOGETHER WITH ELECTROSTATIC TINFOIL ENERGIES
C   ---- EXCLUSION CONTRIBUTION IS CALCULATED IN EXCL_FIELD AND CALLED 
C   ---- ONCE IN EWALDX
C   ---- SELF CONTRIBUTION TO THE FIELD IS 0
C   ---- DIRECT CONTRIBUTION WAS CALCULATED IN FRGFLD
C
C      WRITE(6,*)'FRGFLD1 ',PPEFX(1),PPEFY(1),PPEFZ(1)
C      WRITE(6,*)'FRGFLD2 ',PPEFX(2),PPEFY(2),PPEFZ(2)
C      WRITE(6,*)'FRGFLD3 ',PPEFX(3),PPEFY(3),PPEFZ(3)
C
         CALL RECIPR_FIELD
      END IF
      DO LEFP = 1, NPTTPT
         EFLD(1,LEFP) = PPEFX(LEFP)
         EFLD(2,LEFP) = PPEFY(LEFP)
         EFLD(3,LEFP) = PPEFZ(LEFP)
         IF (IFEWLD .AND. EWLDPL) THEN
            DO I=1,3
               EFLD(I,LEFP) = EFLD(I,LEFP)+RCPFLD(I,LEFP)
     *                      +TNFLD(I)
     *                      +EXCLFLD(I,LEFP)
            END DO
C     
C   ---- DIRTY THINGS
            PPEFX(LEFP) = EFLD(1,LEFP)
            PPEFY(LEFP) = EFLD(2,LEFP)
            PPEFZ(LEFP) = EFLD(3,LEFP)
         END IF        
      ENDDO
C
C      WRITE(6,*)'TOTAL STATIC FIELD ',(EFLD(I,1),I=1,3),
C     *           (EFLD(I,2),I=1,3),(EFLD(I,3),I=1,3)
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)
      IF(IDPUNC.EQ.1) RETURN
C
C     --- CALCULATE POLARIZATION ENERGY ---
C
      IUP=0
      DO I=1,NUP
         IUP = IUP + NPPTS(I)
      ENDDO
      ILOW=1
      DO I=1,NLOW-1
         ILOW = ILOW + NPPTS(I)
      ENDDO
      POLNRG = ZERO
      DO II=ILOW,IUP
         POLNRG = POLNRG + DIND(1,II)*PPEFX(II)
     *                   + DIND(2,II)*PPEFY(II)
     *                   + DIND(3,II)*PPEFZ(II)
      ENDDO
      POLNRG = -0.5D+00*POLNRG
      IF(GOPARR) CALL DDI_GSUMF(677,POLNRG,1)
C
C     --- GRAND INDUCED DIPOLE ---
C
      IF (RUNTYP.EQ.ENERGY .OR. RUNTYP.EQ.GRADIENT) THEN
         TMUX = ZERO
         TMUY = ZERO
         TMUZ = ZERO
         DO LEFP=1, NPTTPT
           TMUX = TMUX + DIND(1,LEFP)
           TMUY = TMUY + DIND(2,LEFP)
           TMUZ = TMUZ + DIND(3,LEFP)
         ENDDO
         TMUXD = TMUX*DEBYE
         TMUYD = TMUY*DEBYE
         TMUZD = TMUZ*DEBYE
         IF(MASWRK)
     *   WRITE(IW,9000) TMUX, TMUXD, TMUY, TMUYD, TMUZ, TMUZD
      ENDIF
C
C     WRITE EFATRM/EFBTRM,EFCHG,EFDIP,AND EFQAD BACK.
C
      IF((IEFP2.EQ.1).AND.(IEFPOL.EQ.1)) THEN
        DO I=1,MXFGPT
C         WRITE MP2 DATA FROM TMP BACK TO COMMON BLOCK
          EFATRM(I)=EFATMP(I)
          EFBTRM(I)=EFBTMP(I)
          EFCHG(1,I)=CHGTMP(1,I)
          EFCHG(2,I)=CHGTMP(2,I)
          EFDIP(1,I)=DIPTMP(1,I)
          EFDIP(2,I)=DIPTMP(2,I)
          EFDIP(3,I)=DIPTMP(3,I)
          EFQAD(1,I)=QADTMP(1,I)
          EFQAD(2,I)=QADTMP(2,I)
          EFQAD(3,I)=QADTMP(3,I)
          EFQAD(4,I)=QADTMP(4,I)
          EFQAD(5,I)=QADTMP(5,I)
          EFQAD(6,I)=QADTMP(6,I)
        ENDDO
      END IF
C
      RETURN
C
 9000 FORMAT(/,17X,'INDUCED DIPOLE',4X,'ATOMIC UNITS',3X,'DEBYE'/
     *       17X,14('-'),4X,12('-'),3X,5('-')/
     *       21X,'X(IND)',9X,F8.5,3X,F8.5/
     *       21X,'Y(IND)',9X,F8.5,3X,F8.5/
     *       21X,'Z(IND)',9X,F8.5,3X,F8.5)
C
      END
C
C*MODULE EFDRVR  *DECK EFPPOL2
      SUBROUTINE EFPPOL2(EFLD,EFADD,POLNRG)
C
C --------------------------------------------------------------
C --- THIS ROUTINE CALCULATES FRAGMENT-FRAGMENT POLARIZATION ---
C ---   ENERGIES IN THE ABSENCE OF AN AB INITO MOLECULE      ---
C --------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXPT=2000, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO=0.0D+00, DEBYE=2.541766D+00)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      CHARACTER*8 DNAME,DPSNAM
C
      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 /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
     *                NPRTGO,IDPUNC,IGOFLG
      COMMON /POL2  / CHGHF(2,MXFGPT),DIPHF(3,MXFGPT),QADHF(6,MXFGPT),
     *                AMULTHF(MXFGPT),BMULTHF(MXFGPT),
     *                AMULT2HF(MXFGPT),BMULT2HF(MXFGPT),
     *                EFPOL2(9,MXFRG),DIND2(3,MXFGPT),DINDD2(3,MXFGPT),
     *                DCORD(3,MXPT),CDIS(3,MXFGPT),
     *                CLPD(4*MXFGPT),ZLPD(4*MXFGPT),
     *                NLPD(4*MXFGPT),KFD(MXFGPT),KLD(MXFGPT),
     *                DNAME(MXPT),DPSNAM(MXFGPT)
      COMMON /INDDP2/ TMUX2,TMUY2,TMUZ2,TMUXD2,TMUYD2,TMUZD2,MINMEM2
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /EFPAR / LOCNFG,NLOW,NUP
C
      DIMENSION EFLD(3,NPTTPT),EFADD(3,NPTTPT)
C
      DATA ENERGY  /8HENERGY  /,GRADIENT/8HGRADIENT/
C
C     IF(MASWRK) THEN
C        WRITE(IW,*) 'EFPPOL2: MP2/EFP1 POLARIZATION ENERGY ROUTINE'
C        WRITE(IW,*) 'THE INPUT EFPOL2 MATRIX IS'
C        CALL PRSQ(EFPOL2,NFRG,3,3)
C     END IF
C
      CALL FRGFLD
      DO LEFP = 1, NPTTPT
         EFLD(1,LEFP) = PPEFX(LEFP)
         EFLD(2,LEFP) = PPEFY(LEFP)
         EFLD(3,LEFP) = PPEFZ(LEFP)
      ENDDO
      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)
      IF(IDPUNC.EQ.1) RETURN
C
C     --- CALCULATE 2ND POLARIZATION ENERGY ---
C
      IUP=0
      DO I=1,NUP
         IUP = IUP + NPPTS(I)
      ENDDO
      ILOW=1
      DO I=1,NLOW-1
         ILOW = ILOW + NPPTS(I)
      ENDDO
      POLNRG = ZERO
      DO II=ILOW,IUP
         POLNRG = POLNRG + DIND2(1,II)*PPEFX(II)
     *                   + DIND2(2,II)*PPEFY(II)
     *                   + DIND2(3,II)*PPEFZ(II)
      ENDDO
      POLNRG = -0.5D+00*POLNRG
      IF(GOPARR) CALL DDI_GSUMF(677,POLNRG,1)
C
C     --- GRAND 2ND INDUCED DIPOLE ---
C
      IF (RUNTYP.EQ.ENERGY .OR. RUNTYP.EQ.GRADIENT) THEN
         TMUX2 = ZERO
         TMUY2 = ZERO
         TMUZ2 = ZERO
         DO LEFP=1, NPTTPT
           TMUX2 = TMUX2 + DIND2(1,LEFP)
           TMUY2 = TMUY2 + DIND2(2,LEFP)
           TMUZ2 = TMUZ2 + DIND2(3,LEFP)
         ENDDO
         TMUXD2 = TMUX2*DEBYE
         TMUYD2 = TMUY2*DEBYE
         TMUZD2 = TMUZ2*DEBYE
         IF(MASWRK)
     *   WRITE(IW,9000) TMUX2,TMUXD2, TMUY2,TMUYD2, TMUZ2,TMUZD2
      END IF
C
      RETURN
 9000 FORMAT(/,17X,'INDUCED DIPOLE',4X,'ATOMIC UNITS',3X,'DEBYE'/
     *       17X,14('-'),4X,12('-'),3X,5('-')/
     *       21X,'X(IND)',9X,F8.5,3X,F8.5/
     *       21X,'Y(IND)',9X,F8.5,3X,F8.5/
     *       21X,'Z(IND)',9X,F8.5,3X,F8.5)
C
      END
C
C*MODULE EFDRVR  *DECK EFGRAD
      SUBROUTINE EFGRAD
C
C ---------------------------------------------------------------
C --- THIS ROUTINE CALCULATES THE GRADIENT ASSOCIATED WITH    ---
C --- A FRAGMENT ONLY RUN, I.E., FRAGMENTS IN THE ABSENCE OF  ---
C ---                  AN AB INITIO MOLECULE                  ---
C ---------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000,
     *           MXPT=2000, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (MAXDEN=25*MXATM, LENABC=2000)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,ISEPS,USEPS
C
      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 /FMCOM / X(1)
      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 /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /ISEPS / ISEPS, USEPS
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      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 /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
     *                QDEN(MAXDEN),AR(LENABC),
     *                NSPA,NPS,NPS2,NDEN,NPSPHER
      COMMON /EFPMP2/ MP2TYP
C
      IF(NGLEVL.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) 'NO NUMERICAL GRADIENT WITH FRAGMENTS'
         CALL ABRT
      END IF
C
      IF(IEFP2.EQ.1 .OR. MP2TYP.EQ.1) THEN
        NFRPTS=NMTTPT+NPTTPT+NRTTPT+NTPATM+NTMO+NDTTPT
      ELSE
        NFRPTS=NMTTPT+NPTTPT+NRTTPT+NTPATM+NTMO
      END IF
C
      CALL VALFM(LOADFM)
      LDRG = LOADFM + 1
      LAST = LDRG   + NAT*NAT
C          EXTRA STORAGE FOR SOME SOLVATION MODELS
                NDIMCD=1
      IF(ISEPS) NDIMCD=NAT+NPS
      LEF3   = LAST
      LCOSDG = LEF3   + MAX(1,3*NFRPTS)
      LWORK  = LCOSDG + NDIMCD*NDIMCD
      LAST   = LWORK  + MAX(1,10*NMTTPT)
                    MXTMP = 1
      IF(NFRG.GT.0) MXTMP = MXFGPT
      LEFATMP = LAST    + 3*MXTMP
      LEFBTMP = LEFATMP + 3*MXTMP
      LCHGTMP = LEFBTMP + 3*MXTMP
      LDIPTMP = LCHGTMP + 3*MXTMP
      LQUATMP = LDIPTMP + 3*MXTMP
      LINDTMP = LQUATMP + 3*MXTMP
      LIDDTMP = LINDTMP + 3*MXTMP
      LEF3FRG = LIDDTMP + 3*MXTMP
      LATORQN = LEF3FRG + 3*MXTMP
      LAST    = LATORQN + MAX(1,3*NFRG)
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      CALL VCLR(DEF,1,3*MXFGPT)
      CALL VCLR(ATORQ,1,3*NFRG)
      CALL VNNDER(X(LDRG),NAT,X(LEF3),NFRPTS,DUMMY,1,X(LWORK),
     *            X(LEFATMP),X(LEFBTMP),X(LCHGTMP),X(LDIPTMP),
     *            X(LQUATMP),X(LINDTMP),X(LIDDTMP),X(LEF3FRG),
     *            X(LATORQN),MXTMP)
      CALL EFFT
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE EFDRVR  *DECK MDX
C      SUBROUTINE MDX
C      WRITE(6,*) 'MDX CALLED BY MISTAKE'
C      CALL ABRT
C      STOP
C      END
C*MODULE EFDRVR  *DECK POTGRAD
      SUBROUTINE POTGRAD
      WRITE(6,*) 'POTGRAD CALLED BY MISTAKE'
      CALL ABRT
      STOP
      END
C*MODULE EFDRVR  *DECK POTNRG
      SUBROUTINE POTNRG(METHOD)
      IF(METHOD.LT.0) RETURN
      WRITE(6,*) 'POTNRG CALLED BY MISTAKE'
      CALL ABRT
      STOP
      END
C*MODULE EFDRVR  *DECK POTMAX
      SUBROUTINE POTMAX(NFRG,DEFT,TORQ,MAXGRD,RMSGRD)
      DOUBLE PRECISION DEFT(3,*),TORQ(3,*),MAXGRD,RMSGRD
      IF(NFRG.LE.0) RETURN
C
      DEFT(1,1) = 0.0D+00
      TORQ(1,1) = 0.0D+00
      MAXGRD = 0.0D+00
      RMSGRD = 0.0D+00
      WRITE(6,*) 'POTMAX CALLED BY MISTAKE'
      CALL ABRT
      STOP
      END
C
C*MODULE EFDRVR  *DECK EFPDYN
      SUBROUTINE EFPDYN(EDISC8,DYNTR,MXPPTS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXFRG=1050, MXIFRQ=12, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      PARAMETER (ZERO=0.0D+00)
      PARAMETER (THREE=3.0D+00, FOUR=4.0D+00, ONE=1.0D+00, TWO=2.0D+00)
      PARAMETER (SIX=6.0D+00)
      PARAMETER (B=1.5D+00)
      PARAMETER (FACT1= 0.72086099022968040154D-02)
      PARAMETER (FACT2= 0.17697067815034886394D-01)
      PARAMETER (FACT3= 0.30660908596251749739D-01)
      PARAMETER (FACT4= 0.48381293256249884995D-01)
      PARAMETER (FACT5= 0.74878830420650517080D-01)
      PARAMETER (FACT6= 0.11806515901361630228D+00)
      PARAMETER (FACT7= 0.19535413832209084204D+00)
      PARAMETER (FACT8= 0.35055692324483221824D+00)
      PARAMETER (FACT9= 0.71577113554429568336D+00)
      PARAMETER (FACT10=1.81409759976323969729D+00)
      PARAMETER (FACT11=6.97923445114870823247D+00)
      PARAMETER (FACT12=83.2480938829658453917D+00)
C
      CHARACTER*8 POLNAM,DPOLNAM
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /EFPAR / LOCNFG,NLOW,NUP
      COMMON /EFPPBC/ XBOX,YBOX,ZBOX,XCUT,YCUT,ZCUT,RCUT
     *                ,SWR1,SWR2,ISWF
      COMMON /EFPGES/ DESFRG(6,MXFRG)
      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 /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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION  DYNTR(MXIFRQ*MXPPTS,NFRG)
C
C     ------------------------------------------------------------------
C        THIS ROUTINE CALCULATES FRAGMENT-FRAGMENT DISPERSION (EFP2)
C        ENERGY AND GRADIENTS IN THE ABSENCE OF AN AB INITO MOLECULE
C     ------------------------------------------------------------------
C
      EDISC8 = ZERO
      EDISD  = ZERO
C
C     -- TRACE OF DYNAMIC ALPHA IN THE FORM ALPHA(NDPPTS,NFRG) ---
C
      ID=0
      DO I=1,NFRG
        DO J=1,NDPPTS(I)
           ID=ID+1
           DYNTR(J,I) = (ONE/THREE)*(EFDPOL(1,ID) + EFDPOL(2,ID) +
     *                   EFDPOL(3,ID))
        ENDDO
      ENDDO
C
      ID1=0
      IF (GOPARR) THEN
         DO I=1,NLOW-1
            ID1 = ID1 + NPPTS(I)
         ENDDO
      ENDIF
C
      DO 280 IM = NLOW,NUP
         INPPTS=NPPTS(IM)
         COMX1 = EFCENT(1,IM)
         COMY1 = EFCENT(2,IM)
         COMZ1 = EFCENT(3,IM)
         DO 270 IN=1,NPPTS(IM)
            ID1=ID1+1
            ALP1=DYNTR(IN,IM)
            ALP2=DYNTR(IN+INPPTS,IM)
            ALP3=DYNTR(IN+2*INPPTS,IM)
            ALP4=DYNTR(IN+3*INPPTS,IM)
            ALP5=DYNTR(IN+4*INPPTS,IM)
            ALP6=DYNTR(IN+5*INPPTS,IM)
            ALP7=DYNTR(IN+6*INPPTS,IM)
            ALP8=DYNTR(IN+7*INPPTS,IM)
            ALP9=DYNTR(IN+8*INPPTS,IM)
            ALP10=DYNTR(IN+9*INPPTS,IM)
            ALP11=DYNTR(IN+10*INPPTS,IM)
            ALP12=DYNTR(IN+11*INPPTS,IM)
            X1 = EFP(1,ID1)
            Y1 = EFP(2,ID1)
            Z1 = EFP(3,ID1)
C
            ID2=0
            DO 260 JM = 1,NFRG
               IF (IM.GE.JM)THEN
                 ID2=ID2+NPPTS(JM)
                 GO TO 260
               END IF
C
               COMX2 = EFCENT(1,JM)
               COMY2 = EFCENT(2,JM)
               COMZ2 = EFCENT(3,JM)
               COMX = COMX1 - COMX2
               COMY = COMY1 - COMY2
               COMZ = COMZ1 - COMZ2
               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
               SWF   = 1.0D+00
               SWFDX = 0.0D+00
               SWFDY = 0.0D+00
               SWFDZ = 0.0D+00
C
               IF(X12.GT.XCUT.OR.Y12.GT.YCUT.OR.Z12.GT.ZCUT)THEN
                 ID2=ID2+NPPTS(JM)
                 GOTO 260
               END IF
               R12  = X12*X12+Y12*Y12+Z12*Z12
               IF(R12.GT.RCUT)THEN
                 ID2=ID2+NPPTS(JM)
                 GOTO 260
               END IF
C
               IF(R12.GT.SWR2*SWR2) THEN
                 ID2=ID2+NPPTS(JM)
                 GOTO 260
               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
                    SWFDX = DSWF*(EFCENT(1,IM)-EFCENT(1,JM)-PBCX)
                    SWFDY = DSWF*(EFCENT(2,IM)-EFCENT(2,JM)-PBCY)
                    SWFDZ = DSWF*(EFCENT(3,IM)-EFCENT(3,JM)-PBCZ)
                 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)
                    DUMY2 =-0.50D+00 * SIN(DUMY) * PIL / R12
                    SWFDX = DUMY2*(EFCENT(1,IM)-EFCENT(1,JM)-PBCX)
                    SWFDY = DUMY2*(EFCENT(2,IM)-EFCENT(2,JM)-PBCY)
                    SWFDZ = DUMY2*(EFCENT(3,IM)-EFCENT(3,JM)-PBCZ)
                 END IF
               END IF
C
               JNPPTS=NPPTS(JM)
               DO 250 JP = 1,NPPTS(JM)
                  ID2=ID2+1
                  ALPN1=DYNTR(JP,JM)
                  ALPN2=DYNTR(JP+JNPPTS,JM)
                  ALPN3=DYNTR(JP+2*JNPPTS,JM)
                  ALPN4=DYNTR(JP+3*JNPPTS,JM)
                  ALPN5=DYNTR(JP+4*JNPPTS,JM)
                  ALPN6=DYNTR(JP+5*JNPPTS,JM)
                  ALPN7=DYNTR(JP+6*JNPPTS,JM)
                  ALPN8=DYNTR(JP+7*JNPPTS,JM)
                  ALPN9=DYNTR(JP+8*JNPPTS,JM)
                  ALPN10=DYNTR(JP+9*JNPPTS,JM)
                  ALPN11=DYNTR(JP+10*JNPPTS,JM)
                  ALPN12=DYNTR(JP+11*JNPPTS,JM)
C
                  X2 = EFP(1,ID2)
                  Y2 = EFP(2,ID2)
                  Z2 = EFP(3,ID2)
                  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)
                  R4 = R2*R2
                  R6 = R4*R2
                  RB = B*R
                  EX1 = EXP(-RB)
                  RB2 = RB*RB
                  RB3 = RB*RB2
                  RB4 = RB2*RB2
                  RB5 = RB*RB4
                  RB6 = RB*RB5
C
                  C6=
     *             + FACT1*ALP1*ALPN1
     *             + FACT2*ALP2*ALPN2
     *             + FACT3*ALP3*ALPN3
     *             + FACT4*ALP4*ALPN4
     *             + FACT5*ALP5*ALPN5
     *             + FACT6*ALP6*ALPN6
     *             + FACT7*ALP7*ALPN7
     *             + FACT8*ALP8*ALPN8
     *             + FACT9*ALP9*ALPN9
     *             + FACT10*ALP10*ALPN10
     *             + FACT11*ALP11*ALPN11
     *             + FACT12*ALP12*ALPN12
C
C                 - F6 IS DAMPING FUNCTION FOR DISPERSION BASED ON
C                   TANG AND TOENNIES FORMULA
C
                  F6 = ONE-(EX1*(ONE + RB + RB2/TWO +
     *                           RB3/SIX + RB4/24.0D+00 +
     *                           RB5/120.0D+00 +RB6/720.0D+00))
C
C                 - EDISD IS DAMPED DISPERSION ENERGY
C
                  EPAIR = - F6*C6/R6
                  EDISD = EDISD + SWF*EPAIR
C
C                 - GRADIENT OF THE DISPERSION ENERGY
C                 DE/DR=FOUR/THREE*(-6*C6*F6/R**7 + DF6/DR * C6/R**6)
C
                  CON1=-(FOUR/THREE)*(C6/R6)*
     *                  ((SIX*F6/R2)-(B*RB6*EX1)/(R*720.0D+00))
C
                  CONX=CON1*X
                  CONY=CON1*Y
                  CONZ=CON1*Z
C
                  EPAIR=EPAIR*FOUR/THREE
                  DESFRG(1,JM) = DESFRG(1,JM) - EPAIR*SWFDX
                  DESFRG(2,JM) = DESFRG(2,JM) - EPAIR*SWFDY
                  DESFRG(3,JM) = DESFRG(3,JM) - EPAIR*SWFDZ
                  DESFRG(1,IM) = DESFRG(1,IM) + EPAIR*SWFDX
                  DESFRG(2,IM) = DESFRG(2,IM) + EPAIR*SWFDY
                  DESFRG(3,IM) = DESFRG(3,IM) + EPAIR*SWFDZ
C
                  DESFRG(1,JM) = DESFRG(1,JM) + SWF*CONX
                  DESFRG(2,JM) = DESFRG(2,JM) + SWF*CONY
                  DESFRG(3,JM) = DESFRG(3,JM) + SWF*CONZ
                  DESFRG(4,JM) = DESFRG(4,JM) + SWF*
     *                      (CONZ*(EFP(2,ID2)-EFCENT(2,JM))
     *                      -CONY*(EFP(3,ID2)-EFCENT(3,JM)))
                  DESFRG(5,JM) = DESFRG(5,JM) + SWF*
     *                      (CONX*(EFP(3,ID2)-EFCENT(3,JM))
     *                      -CONZ*(EFP(1,ID2)-EFCENT(1,JM)))
                  DESFRG(6,JM) = DESFRG(6,JM) + SWF*
     *                      (CONY*(EFP(1,ID2)-EFCENT(1,JM))
     *                      -CONX*(EFP(2,ID2)-EFCENT(2,JM)))
                  DESFRG(1,IM) = DESFRG(1,IM) - SWF*CONX
                  DESFRG(2,IM) = DESFRG(2,IM) - SWF*CONY
                  DESFRG(3,IM) = DESFRG(3,IM) - SWF*CONZ
                  DESFRG(4,IM) = DESFRG(4,IM) - SWF*
     *                      (CONZ*(EFP(2,ID1)-EFCENT(2,IM))
     *                      -CONY*(EFP(3,ID1)-EFCENT(3,IM)))
                  DESFRG(5,IM) = DESFRG(5,IM) - SWF*
     *                      (CONX*(EFP(3,ID1)-EFCENT(3,IM))
     *                      -CONZ*(EFP(1,ID1)-EFCENT(1,IM)))
                  DESFRG(6,IM) = DESFRG(6,IM) - SWF*
     *                      (CONY*(EFP(1,ID1)-EFCENT(1,IM))
     *                      -CONX*(EFP(2,ID1)-EFCENT(2,IM)))
 250           CONTINUE
 260        CONTINUE
 270     CONTINUE
 280  CONTINUE
C
      IF(GOPARR) THEN
         CALL DDI_GSUMF(677,EDISD,1)
      END IF
      EDISC8=EDISD+EDISD/THREE
C
      RETURN
      END
C*MODULE EFDRVR  *DECK POLPCM
      SUBROUTINE POLPCM(EFLD,EFADD,QSE,QSN,
     *                  XCTS,YCTS,ZCTS,AS,
     *                  XE,YE,ZE,RE,ISPHE,LIST,POLNRG)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      CHARACTER*8 FRGNME,POLNAM,DPOLNAM
C
      PARAMETER (MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /EFPPAR/ EFP(3,MXFGPT),EFPOL(9,MXFGPT),
     *                EFDP(3,MXIFRQ*MXFGPT),EFDPOL(9,MXIFRQ*MXFGPT),
     *                ENO,DIND(3,MXFGPT),DINDD(3,MXFGPT),
     *                POLNAM(MXFGPT),DPOLNAM(MXIFRQ*MXFGPT),
     *                POLSCR(MXFGPT),IPLSCR,INDSCR
      COMMON /FMCOM / X(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 /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMEFP/ PB_MUL,PC_MUL,PB_POL,PC_POL,P_FF,P_NF,P_FN,P_I,
     *                P_IBIS,P_J,P_NUCC,P_NUCCBIS,ENPCM
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IKREP,IEF,IP_F,NFMOPCM
      COMMON /PCMPRT/ GCAVP,GCAVS,GDISP,GREP,EHFGAS
      COMMON /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
     *                NPRTGO,IDPUNC,IGOFLG
C
      DIMENSION  EFLD(3,NPTTPT),EFADD(3,NPTTPT),
     *           QSE(*),QSN(*),XCTS(*),YCTS(*),ZCTS(*),
     *           AS(*),XE(*),YE(*),
     *           ZE(*),RE(*),ISPHE(*),LIST(*)
C
      PARAMETER (ZERO = 0.00D+00, ONE = 1.00D+00,
     *           TOANGS=0.52917724924D+00, ANTOAU=1.0D+00/TOANGS)
      PARAMETER (GLIM=1.0D-03)
C
C     -------------------------------------------
C       -- EFP/PCM (NO QM ATOM) CALCULATION --
C       EFP INDUCED DIPOLE AND PCM ASC ARE
C       ITERATED TO SELF-CONSISTENCY
C       HUI LI
C       IOWA STATE UNIVERSITY
C       JULY 25TH, 2005
C     -------------------------------------------
C
      NPRTBK = NPRINT
      NPRINT = 817
      CALL MAKCVM
      CALL PEDRAM
      NPRINT = NPRTBK
C
      MXDII1 = MXDIIS + 1
C
      CALL VALFM(LOADFM)
      LQPOT1 = LOADFM + 1
      LQPOT2 = LQPOT1 + NTS
      LQ0    = LQPOT2 + NTS
      LQ1    = LQ0    + NTS
      LQ2    = LQ1    + NTS
      LQ3    = LQ2    + NTS
      LD0    = LQ3    + NTS
      LQA    = LD0    + NTS
      LDIMAT = LQA    + NTS
      LQREP  = LDIMAT + MXDII1*MXDII1
      LTMP0  = LQREP  + 2*NTS*MXDII1
      LTMP1  = LTMP0  + NTS*3
      LTMP2  = LTMP1  + MXDII1
      LIPVT  = LTMP2  + MXDII1*MXDII1
      LPOTTMP= LIPVT  + MXDII1
      LRMUL  = LPOTTMP+ NTS
      LAST   = LRMUL  + MXSP*10
      NEED   = LAST   - LOADFM -1
      CALL GETFM(NEED)
C
C
C     -- EFP STATIC MULTIPOLE POTENTIAL AT EACH TESSERA
C
      IF(IEFC+IEFD+IEFQ.GT.0) THEN
         DO ITS=1,NTS
            XI=XCTS(ITS)
            YI=YCTS(ITS)
            ZI=ZCTS(ITS)
            QADD=ZERO
            DO I=1,NMTTPT
               ALFA=EFATRM(I)
               BETA=EFBTRM(I)
               XX=EFQAD(1,I)
               YY=EFQAD(2,I)
               ZZ=EFQAD(3,I)
               XY=EFQAD(4,I)
               XZ=EFQAD(5,I)
               YZ=EFQAD(6,I)
               DUM=XX+YY+ZZ
               QXX=(3.0D+00*XX-DUM)/2.0D+00
               QYY=(3.0D+00*YY-DUM)/2.0D+00
               QZZ=(3.0D+00*ZZ-DUM)/2.0D+00
               QXY=1.50D+00*XY
               QXZ=1.50D+00*XZ
               QYZ=1.50D+00*YZ
               DX = XCTS(ITS)-EFC(1,I)
               DY = YCTS(ITS)-EFC(2,I)
               DZ = ZCTS(ITS)-EFC(3,I)
               R2 =  (EFC(1,I)-XI)**2
     *              +(EFC(2,I)-YI)**2
     *              +(EFC(3,I)-ZI)**2
               R  = SQRT(R2)
               FACT = ONE-(BETA*EXP(-ALFA*R2))
               QADD = QADD + FACT*EFCHG(1,I)/R + EFCHG(2,I)/R
     *                     + FACT*
     *                       (EFDIP(1,I)*DX
     *                       +EFDIP(2,I)*DY
     *                       +EFDIP(3,I)*DZ)/(R*R2)
     *                     + FACT*
     *                       (QXX*3.0D+00*DX*DX
     *                       +QYY*3.0D+00*DY*DY
     *                       +QZZ*3.0D+00*DZ*DZ
     *                       +QXY*6.0D+00*DX*DY
     *                       +QXZ*6.0D+00*DX*DZ
     *                       +QYZ*6.0D+00*DY*DZ)
     *                       /(3.0D+00*R2*R2*R)
            ENDDO
            X(LQPOT1+ITS-1)=-QADD
         ENDDO
      END IF
C
C
      CALL REGINI(ISPHE,LIST)
C
C
C     -- FIELD AT POLARIZABILITY POINTS DUE TO
C        MULTIPOLE POINTS IN OTHER FRAGMENTS
C        NOTE: THIS IS DONE ONCE AND REPEATEDLY
C              USED TO INDUCE EFP DIPOLES
C
      CALL FRGFLD
C
C
C     IF THERE IS NO DIPOLE POLARIZABILITY TENSOR,
C     NO ITERATION IS NEEDED.
C
      MXCALL = 1
      IF(IEFP.GT.0) MXCALL = 40
C
      DO 600 ICALL = 1, MXCALL
C        USUALLY 10~15 ITERATIONS ARE REQUIRED.
C
C        -- 1. ADD INDUCED DIPOLE POTENTIAL TO EACH TESSERA
         DO ITS=1,NTS
            XI=XCTS(ITS)
            YI=YCTS(ITS)
            ZI=ZCTS(ITS)
            QADD=ZERO
            DO LEFP=1,NPTTPT
               XL=EFP(1,LEFP)
               YL=EFP(2,LEFP)
               ZL=EFP(3,LEFP)
               X1 = XI - XL
               Y1 = YI - YL
               Z1 = ZI - ZL
               R2=X1*X1 + Y1*Y1 + Z1*Z1
               R = SQRT(R2)
               R3 = R*R2
               QADD= QADD + (DIND(1,LEFP)*X1
     *                    +  DIND(2,LEFP)*Y1
     *                    +  DIND(3,LEFP)*Z1)/R3
            ENDDO
            X(LQPOT2+ITS-1)=X(LQPOT1+ITS-1)-QADD
         ENDDO
C
C        -- 2. COMPUTE PCM SURFACE CHARGE ON EACH TESSERA
C
C        ITERATIVE CYCLE (STEP 1)
         NTSPAR=(NTS-1)/NPROC+1
         CALL ASCCYC(1,ICALL,X(LQPOT2),X(LQ0),X(LQ1),X(LD0),
     *              X(LQA),X(LDIMAT),X(LQREP),X(LTMP0),X(LTMP1),
     *              X(LTMP2),X(LIPVT),X(LRMUL),
     *              MXDII1,NTSPAR,XCTS,YCTS,ZCTS,AS,
     *              XE,YE,ZE,RE,ISPHE,LIST)
C
C        C-PCM ONLY NEEDS STEP 1.
         IF(IEF.EQ.10) THEN
           DO ITS=1,NTS
              QSE(ITS)    =X(LQ1+ITS-1)
              X(LQ0+ITS-1)=X(LQ1+ITS-1)
           ENDDO
           GOTO 500
         END IF
C
C        ITERATIVE CYCLE (STEP 2)
C
C        WRITE THE INTERMEDIATE POTENTIAL, THE STEP 2 INITIAL GUESS
C        AND COPY THE INTERMEDIATE CHARGES
C
         FACIP=-1.0D+00/(EPS-1.0D+00)
         FACIG=-1.0D+00/EPS
         DO ITS=1,NTS
            X(LQ0+ITS-1)=      X(LQ2+ITS-1)
            X(LQ2+ITS-1)=FACIP*X(LQ1+ITS-1)
            X(LQ3+ITS-1)=      X(LQ1+ITS-1)
            IF(ICALL.EQ.1) X(LQ0+ITS-1)=FACIG*X(LQ1+ITS-1)
         ENDDO
C
         CALL ASCCYC(2,ICALL,X(LQ2),X(LQ0),X(LQ1),X(LD0),
     *              X(LQA),X(LDIMAT),X(LQREP),X(LTMP0),X(LTMP1),
     *              X(LTMP2),X(LIPVT),X(LRMUL),
     *              MXDII1,NTSPAR,XCTS,YCTS,ZCTS,AS,
     *              XE,YE,ZE,RE,ISPHE,LIST)
C
C        COPY THE CHARGES
C
         DO ITS=1,NTS
            QSE(ITS)=X(LQ3+ITS-1)+X(LQ1+ITS-1)
         ENDDO
C
C        USE Q1 AND Q3 FOR NEXT GUESSES
C
         DO ITS=1,NTS
            X(LQ0+ITS-1)=X(LQ3+ITS-1)
            X(LQ2+ITS-1)=X(LQ1+ITS-1)
         ENDDO
C
 500     CONTINUE
C
C        -- SCALE C-PCM CHARGE BY (EPS-1)/EPS
         IF(IEF.EQ.10) THEN
           DO ITS = 1, NTS
              QSE(ITS) = QSE(ITS)*(EPS-ONE)/EPS
           END DO
         END IF
C
C        TEST PCM ASC FOR CONVERGENCY
C        NOTE: WE USE QSN TO SAVE THE OLD QSE
C              MAKE SURE TO CLEAR QSN AFTERMATH
C
         SDEV=ZERO
         DO ITS=1,NTS
            SDEV=SDEV+(QSE(ITS)-QSN(ITS))**2
            QSN(ITS) = QSE(ITS)
         ENDDO
         SDEV=SQRT(SDEV/NTS)
C           IF(MASWRK) WRITE(IW,*) ' THRES=',THRES,
C    *      ' SDEV=', SDEV
         IF (SDEV.EQ.ZERO) THEN
C           IF(MASWRK) WRITE(IW,*)
C    *      'POLPCM CONVERGED AFTER ',ICALL,'ITERATIONS'
            CALL VCLR(QSN,1,NTS)
            GO TO 700
         END IF
C
C
C        **********************
C        UPDATE INDUCED DIPOLES
C        **********************
C
         IF(IEFP.GT.0) THEN
C
C           -- TOTAL FIELDS AT POL-TENSORS
C
            IPCOUNT = ME-1
            DO 200 LEFP = 1, NPTTPT
C              -- GO PARALLEL!
               IF(GOPARR) THEN
                  EFLD(1,LEFP) = ZERO
                  EFLD(2,LEFP) = ZERO
                  EFLD(3,LEFP) = ZERO
                  IPCOUNT = IPCOUNT + 1
                  IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 200
               END IF
C              -- SOLVENT APPARENT CHARGE FIELD AT THE POINT --
               PQ_X=ZERO
               PQ_Y=ZERO
               PQ_Z=ZERO
               IF(IEFPOL.EQ.1) THEN
                  DO 330 ITS = 1,NTS
                     RR = (EFP(1,LEFP)-XCTS(ITS))**2 +
     *                    (EFP(2,LEFP)-YCTS(ITS))**2 +
     *                    (EFP(3,LEFP)-ZCTS(ITS))**2
                     IF(RR.LT.GLIM .OR. RR.LT.(REFPOL*ANTOAU))GOTO 330
                     RR3 = RR * SQRT(RR)
                     PQ_X = PQ_X + QSE(ITS)*(EFP(1,LEFP)-XCTS(ITS))/RR3
                     PQ_Y = PQ_Y + QSE(ITS)*(EFP(2,LEFP)-YCTS(ITS))/RR3
                     PQ_Z = PQ_Z + QSE(ITS)*(EFP(3,LEFP)-ZCTS(ITS))/RR3
 330              CONTINUE
               END IF
C              -- TOTAL FIELD: EFP+PCM --
               EFLD(1,LEFP) = PQ_X+PPEFX(LEFP)
               EFLD(2,LEFP) = PQ_Y+PPEFY(LEFP)
               EFLD(3,LEFP) = PQ_Z+PPEFZ(LEFP)
 200        CONTINUE
            IF(GOPARR) THEN
               CALL DDI_GSUMF(667,EFLD, 3*NPTTPT)
            END IF
C
C           -- COMPUTE DIND/DINDD
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
            NEEED    = LAST - LOADFM - 1
            CALL GETFM(NEEED)
            CALL DIPIT(EFLD,EFADD,X(LDINDC),X(LDINDDC),X(LDIND2C),
     *                 X(LDINDD2C),X(LDIND_W),X(LDIND_C),X(LDIND2_W),
     *                 X(LDIND2_C),NPTTPT,X(LDIPNEW),X(LDIPNWD),
     *                 X(LEWLD1),X(LEWLD2),NFRG)
            CALL RETFM(NEEED)
            IF(IDPUNC.EQ.1) RETURN
         END IF
C
 600  CONTINUE
 700  CONTINUE
C
C
C     -- COMPUTE PCM INTERACTION
C
      ENPCM  = ZERO
      QET    = ZERO
      DO ITS = 1,NTS
         QET    = QET+QSE(ITS)
         ENPCM  = ENPCM - X(LQPOT1+ITS-1) * QSE(ITS)
      ENDDO
C
C     -- COMPUTE EFP POLARIZATION ENERGY
C
      POLNRG = ZERO
      DO II=1,NPTTPT
         POLNRG = POLNRG - DIND(1,II)*PPEFX(II)
     *                   - DIND(2,II)*PPEFY(II)
     *                   - DIND(3,II)*PPEFZ(II)
      END DO
      POLNRG = POLNRG/2.0D+00
C
C     SOLVATION OUTPUT FOR PCM MODEL
C
      IF(MASWRK) THEN
        WRITE(IW,*) ' '
        WRITE(IW,*) ' NTS    =', NTS
        WRITE(IW,'(A10,F18.10)') 'PCM ASC=',QET
        WRITE(IW,'(A10,F18.10)') 'GELEC  =',ENPCM/2.0D+00
        IF(ICAV.EQ.1)
     *  WRITE(IW,'(A10,F18.10)') 'GCAVP  =',GCAVP/627.509541D+00
        IF(IDISP.EQ.1)THEN
        WRITE(IW,'(A10,F18.10)') 'GDISP  =',GDISP/627.509541D+00
        WRITE(IW,'(A10,F18.10)') 'GREP   =', GREP/627.509541D+00
        END IF
      END IF
C
      CALL RETFM(NEED)
      RETURN
      END
C
C
C*MODULE EFDRVR  *DECK PREFP
      SUBROUTINE PREFP(LU)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      CHARACTER*6 FRGNAM
      CHARACTER*8 PTNAM
C
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRAGMT/ XCRD(3,MXFRG),YCRD(3,MXFRG),ZCRD(3,MXFRG),
     *                PTNAM(3,MXFRG),FRGNAM(MXFRG)
      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)
C
C   ---- LVS 2007: ROUTINE FOR PRINTING IN MD AND MC
C
      DO 30 I=1,NFRG
       IF (MASWRK) THEN
          WRITE(LU,35) I, FRGNAM(I)
          WRITE(LU,34) EFCENT(1,I),EFCENT(2,I),EFCENT(3,I)
          WRITE(LU,33) DEFT(1,I),DEFT(2,I),DEFT(3,I)
          WRITE(LU,31) TORQ(1,I),TORQ(2,I),TORQ(3,I)
          WRITE(LU,36)
       END IF
 30   CONTINUE
C
 35   FORMAT('FRAGMENT #',I6,'  ',A6)
 31   FORMAT('      TORQ',3F16.10)
 33   FORMAT('     FORCE',3F16.10)
 34   FORMAT('    EFCENT',3F16.10)
 36   FORMAT('    ')
C
      RETURN
      END
C
