C 20 Aug 07 - TN  - include module for electrostatic derivative terms
C
C*MODULE FMOESD  *DECK FMOESDER
      SUBROUTINE FMOESDER(l1,l2,da,layfrg,scffrg,idmrec0,dc,wrk1,
     *                    mapi,mapj,map3,maxl1d,irec00,idam,idad,idat,
     *                    orbxch1,enexch1,kodexch,jodexch,iodexch)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      Parameter (MXATM=2000,MXSH=5000,MXGTOT=20000,MXAO=8192,MXGSH=30,
     *           MXG2=MXGSH*MXGSH,MXRT=100,maxl=5,
     *           half=0.5D+00,zero=0.0D+00)
      logical DIRSCF,FDIFF,dirsav,SCHWRZ,PACK2E,GOPARR,DSKWRK,MASWRK,
     *        out,some,bssedim,orbxch,SAVGOP,esppar,nxt,espap,doddcor,
     *        esdder,largepri,LCFLAG,LRINT,LCFLAGs,LRINTs
      LOGICAL ESDER
      logical orbxch1,enexch1,kodexch,jodexch,iodexch
C
      dimension da(*),layfrg(*),scffrg(*),idmrec0(*),karten(0:maxl-1)
      dimension dc(*),wrk1(*),mapi(*),mapj(*),map3(maxl1d,3)
C
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FMCOM / X(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,inttyp,igrdtyp
c     COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),ch(mxgtot),ci(mxgtot),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /SHLEXC/ NORGSH(3),NORGSP(3),IEXCH,NANGM,NGTH(4)
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      common /fmopnt/ lichfg,lmulfg,lidmrec,lfrgnam,llayfrg,lindat,
     *                lnCBS,lfmozan,lfmoc,lfmomas,lizbas,liaglob,libdgh,
     *                liabdfg,ljabdfg,lnCAO,lidxCAO,liaprjo,ljaprjo,
     *                lCoreAO,lOccCor,lshiftb,liodfmo,lfmoda,lfmodb,
     *                lfmoespa,lfmoespb,llocfmo,lscffrg,lrij,lpopmul,
     *                lpopmat,lialoc,lindbd,liatfrg,lindfrg,lindgfrg,
     *                lnatfrg,lnat0frg,lianfrg,lzanfrg,lcfrg,llibish,
     *                llibnsh,llibng,lindatg,lfmobuf(3),lfmode,lnumfrg,
     *                lloctat,liaoglob,lloadm,lfmoge,ldgrid,liodcfmo,
     *                ljob2grp,lfmopg,lemocdr,luntxyz,luntrot,lstonep,
     *                lmapsu,lfrgmul
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
C
      COMMON /GRAD  / DE(3,MXATM)
C
      data karten/1,4,6,10,15/ 
      data rhf/8HRHF     /,RMC/8HMCSCF   /,
     *     dbgfmo/8HDBGFMO  /,dbgme/8HFMOESP  /,debug/8HDEBUG   /
c
c     this routine computes the electrostatic potential (ESP).
c     one electron part is computed in ONEEI. here only two-electron
c     contribution is added.
c     parstat: GroupFull/GroupNone 
c
c     For nbsse=3 monomers and dimers during BSSE runs are in vaccuum
      if(nbsse.eq.3.and.ifmostp.eq.5. or. ifmostp.eq.1) return
C
C     this routine must be called either in monoscf, edimer or etrimer
C
      if (ifmostp.eq.6) then
        write (iw,*) 'ifmostp=6 is not allowed in fmoesder'
        call abrt
      end if
      if (ifmostp.ne.2.and.ifmostp.ne.4.and.ifmostp.ne.8) return
c     if (maswrk) WRITE (6,*) 'NOW INSIDE OF FMOESDER'
CZ
c
      some=(exetyp.eq.debug.or.exetyp.eq.dbgfmo.or.exetyp.eq.dbgme).and.
     *     maswrk
      orbxch=mod(modorb,2).ne.0
c     enexch=mod(modorb/2,2).ne.0
      esppar=mod(modpar/2,2).ne.0.and.goparr
c     .true.  divide fragments
c     .false. divide shells
      lwrkden=lfmobuf(1)
      lwrkesp=lfmobuf(2)
      lwrkesp2=lfmodb
c     for monomer SCF there is no need to adjust irec0 below (and MPPROP may 
c     not have been read anyway). 
      ilay=icurlay
      ifg=icurfg
      if(ifg.eq.0) return
c
      lcflags=lcflag
      lrints=lrint
      lcflag=.false.
      lrint=.false.
c
c     ESP is not needed. This happens when monomer/dimer SCF iterations run.
c     (ESP is computed before such iterations along with 1-el integrals)
c
      CALL DERCHK(NDER)
CZ
      esder = nder.gt.0.and.iand(modgrd,2).ne.0
      if (.not.esder) then
        write(iw,*) 'must debug for esder'
        call abrt
      end if
CZ    esdder=ifmostp.eq.6.and.resdim.ne.0.and.iand(ixesp,32).eq.0.and.
CZ   *       nder.gt.0
      esdder=.true.
CZ
CZ    (1). Density derivatives of ESPs
CZ        Note that OUTPUT DA is DI, delta DIJ or delta DIJK
CZ
c     if (maswrk) write(6,*) 'Before ESDNSDER'
      CALL ESDNSDER(L1,L2,DA,X(LFMODB),DC,X(LWRKDEN),WRK1,
     *     layfrg,idmrec0,X(LFMOESPA),
     *     X(LWRKESP),X(LNUMFRG),mapi,mapj,map3,maxl1d,irec00,
     *     orbxch1,enexch1,kodexch,jodexch,iodexch)
CZ
cz    if (maswrk) then
cz      write (6,*) 'after Density matrix'
cz      do i = 1, l1
cz      kk = i*(i-1)/2
cz      write (6,'(6F12.8)') (da(kk+j),j=1,i)
cz      end do
cz    end if
CZ
      jfg=jcurfg
      lfg=kcurfg
CZ
Cz    if (maswrk) WRITE (6,'(A,6I4)') 'ifg jfg lfg',ifg,jfg,lfg
CZ
      kfg0=1
      kfg1=nfg
      iz=1
CZ    kfact=1
CZ
      respapi=respap(iz)
      resppci=resppc(iz)
      bssedim=nbsse.eq.2.and.ifmostp.eq.5
      if(bssedim) jfg=0
      doddcor=iand(ixesp,1024).ne.0
c
c     SCHWRZ = ISCHWZ.GT.0
      SCHWRZ = .true. 
c     Do not think that you can set SCHWRZ to .false. and get away with it.
c     SCHWRZ also forces integral initialisation that will otherwise not be 
c     done.
c 
      dirsav=dirscf
      dirscf=.true.
c     scftyp1=scftyp
c     if(scftyp.eq.rmc) scftyp1=rhf
      scftyp1=rhf
      out=SCHWRZ.and.maswrk.and.iand(nprfmo,3).eq.0
c
      call vclr(x(lfmoespa),1,l2)
cb    call vclr(x(lfmoespb),1,l2)
c
c
c     save the pristine monomer(dimer) configuration
c
      call monbsr(nat,ich,mul,num,nqmt,ne,na,nb,nshell,ngau,enucr,
     *            nat0,ich0,mul0,num0,nqmt0,ne0,na0,nb0,
     *            ncursh,ngau0,enucr0)
c
c     fix the number of electrons to avoid double counting:
c     exclude charge as otherwise the charge of fragment I will be added
c     to NE first here and then again in makmol below.
c     At present ne is actually not important because numfrg has na.
      ne0c=ne0+ich0
c
      nprsav=nprfmo
      if(nprint.eq.-5.and.iand(nprfmo,3).eq.0) nprfmo=nprfmo+1
c     that is, enforce reduced output if nprint.eq.-5
c     For ESPs the only meaningful usage of loadm is with true esppar.
c     loadm should only be used if full range of fragments is treated
c     (that is, excluding separated dimers).
      loadhf=mod(modpar,2)
      largepri=loadhf.eq.1.and.esppar.and.kfg0.eq.1.and.kfg1.eq.nfg
c     write(6,*) 'wwwlb',largepri
CZ
CZ    (2). BASIS FUNCTION DERIVATIVES OF 1E INTEGRALS AND HELMAN-FEYNMAN TERMS
CZ    DA is obtained in ESDNSDER 
CZ    The contributions are directly stored in FMODE().
CZ
      CALL ES1DER(da,x(lwrkden),x(lwrkesp),L2,0,idam,idad,idat)
cz
      natf1 = natfmo*3
      natf2 = natf1 *2
CZ
c
C     ----- INITIALIZATION FOR PARALLEL WORK -----
C
      NXT = IBTYP.EQ.1
      NEXT  = -1
      kount = -1
c
c     ESP will be computed in the direct fashion
c
c     do 100 kfg=1,nfg
      do 100 kkfg=kfg0,kfg1
        if(largepri) then
          kfg=ixftch(x(lloadm),kkfg)
        else
          kfg=kkfg
        endif
        if(ifg.eq.kfg.or.jfg.eq.kfg.or.lfg.eq.kfg) goto 100
c
        if(resppci.ne.zero.or.respapi.ne.zero.or.resdim.ne.zero) then
          rk=fmodist(ifg,jfg,lfg,kfg)
c         the distance will be refined later for smartr. 
c         write(6,*) 'Distance is',kfg,rk,resppci
CZCZ      if(iz.eq.1.and.rk.gt.resppci.and.resppci.ne.zero) goto 100
          if(rk.gt.resppci.and.resppci.ne.zero) goto 100
CZ
c
         if(ifmostp.eq.2.and.iand(ixesp,4096).ne.0.and.rk.ne.0) goto 100
c
          espap=rk.gt.respapi.and.respapi.ne.zero
        else
          espap=.false.
        endif
c
CZ
CZ      WHEN IS ESPPAR TRUE?????? 03/16/07
CZ      I HAVE TO ASK DMITRI 
CZ
        if(esppar) then
          kount=kount+1
          if(goparr) then
            IF(NXT) THEN
              IF(kount.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
              if(NEXT.ne.kount) goto 100
            else
              if(MOD(kount,NPROC).NE.me) goto 100
            endif
          endif
          savgop=goparr
          goparr=.false.
c         goparr is set to .false. to prevent shell-based work division
c         inside of EXCHNG and FMOESP. 
c         Note that shell-based work division is used for esppar=.false. 
        endif
        klay=min(ilay,layfrg(kfg))
c       bssedim=nbsse.eq.2.and.ifmostp.eq.5.and.kfg.eq.jcurfg
        bssedim=nbsse.eq.2.and.ifmostp.eq.5
        if(some) write(6,*) 'Computing ESD due to frg',kfg,' layer',klay
c
c       add monomer kfg and read its density (now only alpha)
c
        kfgx = kfg
        call makemol(kfgx,0,0,klay,0,nat0,ncursh,ngau0,ne0c,ich0,mul0,
     *               .false.)
        l1k=num-num0
        l2k=(l1k*l1k+l1k)/2
        l3k=l1k*l1k 
c       write(6,*) 'Computing ESP due to frg',kfg,' layer',klay
c      write(6,*) (x(lpopmul+((icurpop-1)*nfg+kfg-1)*maxl1+i-1),i=1,l1k)
c       write(6,*) (x(lpopmat+((icurpop-1)*nfg+kfg-1)*maxnat+i-1),i=1,
c    *              nat-nat0)
c
c       for separated dimers call 1e integrals and quit. 
c       luckily for separated dimers the distance given by fmodist and fmosdist
c       are exactly same (monomer-monomer distance).
CZCZ
cz      IF (LESDIM) THEN
cz        IFLAG = 1
cz        call vclr(x(lfmoespb),1,l2)
CZ
c         set nat1e to jfg atoms
cz        nat1e=nat-nat0 
c         basis set must be reset to ifg monomer to do 1e integrals. 
c         save dimer info
cz        nshs=nshell
cz        nums=num
cz        nats=nat
cz        num=num0
cz        nshell=ncursh
cz        nat=nat0
c         icurfg=0 prevents unwanted recursive calls to FMOESP from ONEEI.
c         icurfg is stored in ncursh, that is unused in oneei.
cz        ncurs=ncursh
cz        ncursh=icurfg 
cz        icurfg=0
CZ
cz        RESPPCS = RESPPC(2)
cz        ! not to calculate approximation to 2e terms
cz        RESPPC(2) = ZERO
cz        IFMOS   = IFMOSTP
cz        IFMOSTP = 6
CZ
c
cz        call oneei
c
CZ
cz        IFMOSTP = IFMOS
cz        RESPPC(2) = RESPPCS
CZ
cz        ncursh=ncurs
cz        nat1e=natfmo
cz        num=nums
cz        nshell=nshs
cz        nat=nats
c         resppci uses point charges - no 2e integrals; factk=half skips
c         upper triangle of ifg,jfg contributions as they are symmetric 
c         add interfragment nuclear repulsion, dEIJ=(EIJ-EI-EJ)/2
cz        n0j=nat0+1
c
c         The condition below ifg.gt.jfg is placed to avoid double counting.
c         FMOESP is called twice, so there would have been two equal 
c         Save 1e potential 
cz        CALL DAXPY(L2,ONE,X(LFMOESPB),1,X(LFMOESPA),1)
cz        if(rk.gt.resppci.and.resppci.ne.zero) goto 100
cz      endif
c

        irec0=idmrec0(klay)
c       hack it up to point to RHF densities for runs with MP2 densities. 
        if(doddcor) then
          irec0rhf=1
          if(irec0.eq.1) irec0rhf=nfg*2+1
          irec0=irec0rhf
        endif
        idmrec0k=kfg+irec0
        if(.not.espap) then
        if(orbxch) then 
c         nak=na-na0
          nak=ishft(ixftch(x(lnumfrg),kfg),-16) 
          if(scffrg(kfg).eq.rmc) then
CZ          call rareads(IDAFMO,x(liodfmo),x(lfmoda),l2k+l3k,idmrec0k,0)
           call rareads(IDAFMO,x(liodfmo),x(lwrkden),l2k+l3k,idmrec0k,0)
          else
CZ          call rareads(IDAFMO,x(liodfmo),x(lfmoda+l2k),l3k,idmrec0k,0)
CZ          call DMTX2(x(lfmoda),x(lfmoda+l2k),nak,l1k,l1k) 
           call rareads(IDAFMO,x(liodfmo),x(lwrkden+l2k),l3k,idmrec0k,0)
            call DMTX2(x(lwrkden),x(lwrkden+l2k),nak,l1k,l1k) 
          endif
        else
CZ        call rareads(IDAFMO,x(liodfmo),x(lfmoda),l2k,idmrec0k,0)
          call rareads(IDAFMO,x(liodfmo),x(lwrkden),l2k,idmrec0k,0)
        endif
c       goto 100
c       write(6,*) 'original Density',kfg,l1k,nak
c       call prtri(X(lfmoda),l1k)
c       CALL SEQOPN(38,'HESSIAN','UNKNOWN',.FALSE.,'UNFORMATTED')
c       write(38) (X(lfmoda+i-1),i=1,l2k)
c       CALL SEQCLO(38,'KEEP')
c       read beta density
        endif
c
c       array sizes are different for each fragment. 
c
        NSH2 = (NSHELL*NSHELL+NSHELL)/2
        CALL BASCHK(LMAX)
        if(lmax.ge.5) then
           if(maswrk) write(iw,*) 'fmo has not been reviewed for h,i'
           call abrt
        end if
c       For some reason GAMESS likes to handle at least L-shells
        NANGM=karten(max(lmax,1))
        MAXG = NANGM**4
        CALL VALFM(LOADFM)
C
        LXINTS= LOADFM + 1
        LGHOND= LXINTS + NSH2
        LDSH  = LGHOND + MAXG
        LDSHb = LDSH   + NSH2
        LDDIJ = LDSHb  + NSH2
        LAST  = LDDIJ  + 49*MXG2
        NEED  = LAST- LOADFM -1
        CALL GETFM(NEED)
c
        ldenp=lpopmul+((icurpop-1)*nfg+kfg-1)*maxl1
        if(espap) then
          ldena=ldenp
          ldenb=ldena
        else
CZCZ      ldena=lfmoda 
          ldena=lwrkden 
          ldenb=lfmodb 
        endif 
c
c       CALL ERIPRE
c       There is no need to call ERIPRE here. FMO code always calls EXCHNG
c       (because SCHWRZ is forced true), ahd this is where initialisation
c       for each (n+1)-mer is done. Note that this initialisation is 
c       important now: it sets arrays for screening and some other things. 
c
c       create density matrix for shell indices
c
        IF(SCHWRZ) THEN
          DUMMY = 0.0D+00
cnb       for gradients, what should be SCFTYP??
c         l1 and l2k seem to be the right choice:
c         l1 is used to skip L1 AOs to get to the external monomer shell
c         l2k gives the size of GVB density matrices for the external monomer
c         write(6,*) 'original Density',kfg,l1k
c         call prtri(X(lfmoda),l1k)
          if(espap) then
c           write(6,*) 'mullikens'
c           call prsq(X(ldena),l1k,1,1)
           CALL SHLDEND(scftyp1,X(ldena),x(ldenb),DUMMY,X(LDSH),L1,1)
c          write(6,*) 'contracted l-Density',kfg,nshell-ncursh
c           call prsq(X(LDSH),nshell-ncursh,1,1)
          else
            CALL SHLDEN(scftyp1,X(ldena),x(ldenb),DUMMY,X(LDSH),IA,
     *                  L1,L2k,NSH2,1)
            if(esdder) then
c            remove all evidence of fragment 2 being involved. 
c            Keep only fragment 1 information and obtain its shell density. 
c            L1 argument is unused below.
             nshsav=nshell
             nshell=ncursh
             ncursh=0
CZCZ        CALL SHLDEN(scftyp1,X(lwrkden),x(ldenb),DUMMY,X(LDSHb),IA,
CZCZ *                  L1,L2k,NSH2,1)
            CALL SHLDEN(scftyp1,da,x(ldenb),DUMMY,X(LDSHb),IA,
     *                  L1,L2k,NSH2,1)
             ncursh=nshell
             nshell=nshsav
            endif
c         write(6,*) 'contracted Density',kfg,nshell-ncursh
c         call prtri(X(LDSH),nshell-ncursh)
          endif
        END IF
c
c       initialise two-electron integrals
c       Pople integrals use a separate coordinate common block initialised
c       above, etc. ?st must be reset so that no phoney restart is done.
c
        ist=1
        jst=1
        kst=1
        lst=1
        call jandk
C
C     ----- EXCHANGE INTEGRALS FOR DIRECT SCF THRESHOLD TESTS -----
c       Computed for the combined system, no need to adjust indices.
C
        NINT=0
        NSCHWZ=0
cnb5    EXCHNG can be adjusted to take advantage of the reduced basis with BSSE.
c       goto 100
c       do iohno=1,1000
        IF(SCHWRZ) CALL EXCHNG(X(LXINTS),X(LGHOND),X(LDDIJ),
     *                         NSH2,MAXG,inttyp)
c       enddo
c       compute 2-el integrals corresponding to placing the original 
c       monomer(dimer) into the Coulomb field of monomer kfg.
c 
c       goto 100
CZZ     if(esdder) call vclr(x(lwrkesp2),1,l2k)
        call vclr(x(lwrkesp2),1,l2k)
CZ
c       call timit(1)
c       do iohno=1,100
CZ
CZ      what I want is not x(lfmoespa), but x(lwrkesp2).
CZ      I will modify fmoesp2 to do so 
CZ
        CALL fmo2ei(SCHWRZ,NINT,NSCHWZ,NSCHWZB,NSCHWNZB,L1,L2,X(LXINTS),
     *              NSH2,X(LGHOND),MAXG,IA,ldena,ldenp,X(lfmoespa),
     *              da,x(lwrkesp2),X(LDSH),X(LDSHb),x(liaglob),
     *              x(lindat),x(lindatg),bssedim,respapi,resppci,espap,
     *              esdder,esder,ifg,jfg,lfg,kfg)
CZCZ *              esdder,ifg,jfg,lfg,kfg)
c       enddo
c       call timit(1)
        if(out) then
          if(NSCHWZB.eq.0) then
            write(iw,9000) kfg,NINT,NSCHWZ
          else
            write(iw,9005) kfg,NINT,NSCHWZ,NSCHWZB
          endif
        endif
c         write(iw,*) 'Non-zero superblocks',NSCHWNZB,NSH2
c
c       write(6,*) 'Density in ESP is'
c       call prtril(x(lfmoda),num0)
CZ
CZ      added by nagata 3/27/07
CZ
CZCZCZ  IF(SCHWRZ) CALL DAWRIT(IDAF,IODA,X(LXINTS),NSH2,54,0)
        IF(SCHWRZ) lfmobuf(3) = LXINTS
CZ
CZCZCZ  CALL RETFM(NEED)
c
        if(esppar) goparr=savgop
CZCZ    if(esdder) then
          CALL DSCAL(L2k,half,X(lwrkesp2),1)
          II=lwrkesp2-1
          DO I=1,L1k
            II = II+I
            X(II) = X(II) + X(II)
          enddo
          if(goparr) call ddi_gsumf(2418,x(lwrkesp2),l2k)
c         write(6,*) 'Having V1',l1k
c         call prtril(X(lwrkesp),l1k)
c         write(6,*) 'Having V2',l1k
c         call prtril(X(lwrkesp2),l1k)
CZCZ      CALL daxpy(L2k,one,X(lwrkesp2),1,X(lwrkesp),1)
CZCZ    endif

CZ      test
c       write(6,*) 'WRKESP2',L1
c       CALL PRTRIL(X(LWRKESP2),L1K)
CZ
c
CZ
CZ      calculate Lagrangian-like term, X
CZ
        ! add square space after triangle space for lfmoda
        CALL VCLR(X(LWRKESP),1,L2K)
CZ
        CALL CPYTSQ(X(LWRKDEN),X(LWRKDEN+L2K),L1K,1)
        CALL TFTRI(X(LWRKESP),X(LWRKESP2),X(LWRKDEN+L2K),
     *             X(LWRKESP+L2K),L1K,L1K,L1K)
        CALL DSCAL(L2K,-HALF,X(LWRKESP),1)
        CALL VCLR(DE,1,NAT*3)
        CALL EXTSDER(X(LWRKESP),X(LWRKESP+L2K),L1K,L2K,.FALSE.,
     *               .FALSE.,L1)
cz      if (maswrk) then
cz        WRITE(6,*) 'test after lagrangian-like term',kfg,nat
cz        do ij = 1, nat
cz          write(6,'(3F12.8)') DE(1,ij),DE(2,ij),DE(3,ij)
cz        end do
cz      end if
CZCZ    CALL PRTRIL(X(LWRKESP),L1K)
CZ
CZ
CZ      CALCULATIONS OF DERIVATIVE TERMS
CZ 
CZ      ELECTRON REPULSION TERM
CZ
CZ      DA is DI, delta DIJ, or delta DIJK 
CZ      X(LWRKDEN)  is the density for external monomer 
CZ      Note that derivatives are stored in DE().
CZ
        CALL ESD2DER(DA,X(LWRKDEN),L1)
        IF(GOPARR) CALL DDI_GSUMF(2418,DE,NAT*3)
cz      if (maswrk) then
cz        WRITE(6,*) 'test after ESD2DER',kfg,nat
cz        do ij = 1, nat
cz          write(6,'(3F12.8)') DE(1,ij),DE(2,ij),DE(3,ij)
cz        end do
cz      end if
        CALL RETFM(NEED)
CZ
        if (maswrk) then
          if (ifmostp.eq.2) then
            call fmodeg(1,x(lfmode),x(lfmopg),x(liaglob))
          elseif (ifmostp.eq.4) then
            call fmodeg(1,x(lfmode+natf1),x(lfmopg),x(liaglob))
            n11 = idam + idad
            n11 = 0 
            if (nbody.gt.2) call fmodeg
     *         (n11,x(lfmode+natf2),x(lfmopg),x(liaglob))
            ! must change
          elseif (ifmostp.eq.8) then
            call fmodeg(1,x(lfmode+natf2),x(lfmopg),x(liaglob))
          endif
        endif
  100 continue
      IF(esppar.and.goparr.and.nxt) CALL DDI_DLBRESET
c
c     restore the pristine monomer(dimer) configuration
c
      call monbsr(nat0,ich0,mul0,num0,nqmt0,ne0,na0,nb0,ncursh,ngau0,
     *           enucr0,nat,ich,mul,num,nqmt,ne,na,nb,nshell,ngau,enucr)
c
c     it is important to reset ncursh so that 2e integrals within the monomer
c     (dimer) during the following SCF are computed properly.
      ncursh=0
c
c     compute 2-el integrals for the original monomer(dimer)
c
      if(ifmostp.ne.6) then
        ist=1
        jst=1
        kst=1
        lst=1
        call jandk
      endif
      nprfmo=nprsav
      dirscf=dirsav
      lcflag=lcflags
      lrint=lrints
c     IF (MPLEVL.EQ.0) THEN
c       CALL DAREAD(IDAF,IODA,DA,L2,16,0)
c     ELSE
c       CALL DAREAD(IDAF,IODA,DA,L2,308,0)
c     END IF
c     if (maswrk) write(6,*) 'Exit of FMOESD',nqmt
      return
 9000 format(/1x,'ESD lfg=',I5,': NZ',I12,' skipped',I10,' blocks.')
 9005 format(/1x,'ESD lfg=',I5,': NZ',I12,' skipped',I10,' blocks',I8,
     *           ' superblocks.')
      end
C*MODULE FMOESD  *DECK ES1DER
      SUBROUTINE ES1DER(DENAB,DTINT,DVINT,L2,MODUS,idam,idad,idat)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL OUT,DBG,NORM
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT,SVDSKW,DEBUGOPT
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192)
      PARAMETER (MXCHRM=1)
C
      DIMENSION DENAB(L2),DTINT(L2,3,NAT),DVINT(L2,3,NAT)
      DIMENSION DIJ(225), IJX(35), IJY(35), IJZ(35),
     *          XV(6,5,5), YV(6,5,5), ZV(6,5,5),
     *          DXV(5,5,5), DYV(5,5,5), DZV(5,5,5)
CZ   *          XS(6,7), YS(6,7), ZS(6,7), XT(6,5), YT(6,5), ZT(6,5),
CZ   *          DXS(5,5), DYS(5,5), DZS(5,5),
CZ   *          DXT(5,5), DYT(5,5), DZT(5,5),
C
      COMMON /CHMGMS/ XCHM(MXCHRM),YCHM(MXCHRM),ZCHM(MXCHRM),
     *                DXELMM(MXCHRM),DYELMM(MXCHRM),DZELMM(MXCHRM),
     *                QCHM(MXCHRM),NCHMAT,KCHRMM
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /FMCOM / X(1)
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      common /fmopnt/ lichfg,lmulfg,lidmrec,lfrgnam,llayfrg,lindat,
     *                lnCBS,lfmozan,lfmoc,lfmomas,lizbas,liaglob,libdgh,
     *                liabdfg,ljabdfg,lnCAO,lidxCAO,liaprjo,ljaprjo,
     *                lCoreAO,lOccCor,lshiftb,liodfmo,lfmoda,lfmodb,
     *                lfmoespa,lfmoespb,llocfmo,lscffrg,lrij,lpopmul,
     *                lpopmat,lialoc,lindbd,liatfrg,lindfrg,lindgfrg,
     *                lnatfrg,lnat0frg,lianfrg,lzanfrg,lcfrg,llibish,
     *                llibnsh,llibng,lindatg,lfmobuf(3),lfmode,lnumfrg,
     *                lloctat,liaoglob,lloadm,lfmoge,ldgrid,liodcfmo,
     *                ljob2grp,lfmopg,lemocdr,luntxyz,luntrot,lstonep,
     *                lmapsu,lfrgmul
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
CZCZ  COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
CZCZ *                MPLEVL,MPCTYP
C
CZ    PARAMETER (MAXDEN=25*MXATM, LENABC=2000, NPPA=1082)
CZ    LOGICAL ISEPS, USEPS
CZ    COMMON /ISEPS / ISEPS, USEPS
CZ    COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
CZ   *                QDEN(MAXDEN),AR(LENABC),
CZ   *                NSPA,NPS,NPS2,NDEN,NPSPHER
CZ    COMMON /SOLVI / IATSP(LENABC+1),N0(2),NP1,NP2,ISKIP
CZ    COMMON /COSMO3/ COSZAN(NPPA),CORZAN(3,NPPA)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00,TWO=2.0D+00)
      PARAMETER (RLN10=2.30258D+00)
      PARAMETER (PI212=1.1283791670955D+00)
      PARAMETER (SQRT3=1.73205080756888D+00)
      PARAMETER (SQRT5=2.23606797749979D+00,SQRT7=2.64575131106459D+00)
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HTVDER   /, GRD1/8HGRD1    /,
     *     RESC/8HRESC    /,ANESC/8HNESC    /
C
      DATA IJX / 1, 2, 1, 1, 3, 1, 1, 2, 2, 1,
     *           4, 1, 1, 3, 3, 2, 1, 2, 1, 2,
     *           5, 1, 1, 4, 4, 2, 1, 2, 1, 3,
     *           3, 1, 3, 2, 2/
      DATA IJY / 1, 1, 2, 1, 1, 3, 1, 2, 1, 2,
     *           1, 4, 1, 2, 1, 3, 3, 1, 2, 2,
     *           1, 5, 1, 2, 1, 4, 4, 1, 2, 3,
     *           1, 3, 2, 3, 2/
      DATA IJZ / 1, 1, 1, 2, 1, 1, 3, 1, 2, 2,
     *           1, 1, 4, 1, 2, 1, 2, 3, 3, 2,
     *           1, 1, 5, 1, 2, 1, 2, 4, 4, 1,
     *           3, 3, 2, 2, 3/
C
C     ----- BASIS FUNCTION DERIVATIVE CONTRIBUTIONS TO GRADIENT -----
C     INTEGRALS ARE OF TYPE <II'/H/JJ> = <II'/T+V/JJ>
C     RESC RUNS STORE INTEGRALS WITHOUT CONTRACTING THEM WITH DENSITY
C
C     MODUS - USED ONLY WITH NESC METHOD AND ONLY FOR KIN. ENERGY T
C           = 0 BULK T GRADIENT
C           = 1 CORRECTION TO S DERIVATIVE (PROPORTIONAL TO T)
C           = 2 CORRECTION TO T DERIVATIVE
C
      DBG = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME)  .AND.  MASWRK
      OUT = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME  .OR.
     *       EXETYP.EQ.GRD1   .OR.  NPRINT.EQ.-3)
      IF(DBG) WRITE(IW,9000)
      IAZ=0
C
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      IF(RMETHOD.EQ.ANESC) THEN
         FSC=ONE/CLIG
         FSC2=FSC*FSC/TWO
      ENDIF
      IF(RMETHOD.EQ.RESC) THEN
         CALL VCLR(DTINT,1,L2*3*NAT)
         CALL VCLR(DVINT,1,L2*3*NAT)
      ENDIF
CZ
CZ    determine the coefficients for FMO n-mer grad 
CZ
      idamm = 0
      idadd = 0
      idatt = 0
      if (ifmostp.eq.2) then
        idamm = idam
        idamm = 1
      elseif (ifmostp.eq.4) then
        if (nbody.gt.1) idadd = 1
        if (nbody.gt.2) idatt = idad + idat ! must change
        if (nbody.gt.2) idatt = 0           ! must change
      elseif (ifmostp.eq.8) then
        idatt = 1 
      else
        write (6,*) 'ifmostp must be 2, 4, or 8'
        call abrt
      endif
CZCZ  if (idamm.eq.0.and.idadd.eq.0.and.idatt.eq.0) return
CZ    write(6,*) 'coef for frag ener in es1der',idamm,idadd,idatt
CZ
C
C INITIALIZE PARALLEL
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT = -1
      MINE = -1
C
C     -D*ELMM- CONTAINS CONTRIBUTION TO THE FORCES ON THE
C     EXTERNAL CHARGES WHICH ARE DUE TO QM MOLECULE.
C
      IF(NCHMAT.NE.0) THEN
         DO 100 ICHARM=1,NCHMAT
            DXELMM(ICHARM)=ZERO
            DYELMM(ICHARM)=ZERO
            DZELMM(ICHARM)=ZERO
  100    CONTINUE
      END IF
CZ    ESD=NFG.NE.0.AND. IFMOSTP.EQ.2.AND.MPLEVL.EQ.0.AND.
CZ   *                  RESDIM.NE.0.AND.IAND(IXESP,32).EQ.0
CZ    NEXTAT=NCHMAT
CZ    IF(ESD) NEXTAT=NATFMO
CZ    NEXTAT=NATFMO
CZ
CZ    set up a parameter in deepest loop    
CZ
      MAXIC = NATFMO
CZ
C
C     FIRST, LOOP OVER THE GENUINE QM MOLECULE (ICHARM=0 PASS),
C     AND THEN ANY EXTERNAL CHARGES WHICH CHARMM MAY BE USING.
C     THE LATTER MAKE A CONTRIBUTION TO THE GRADIENT OF THE ATOMS
C     IN THE QM MOLECULE, AS WELL AS TO THEMSELVES.
C     FMO EXTERNAL ATOMS FROM SEPARATED DIMERS ARE HANDLED SIMILARLY TO
C     CHARMM.
C
CZ
CZ    THIS LOOP IS MOVED 
CZ
CZCZ  DO 2000 ICHARM=0,NEXTAT
CZCZ    IF(ESD.AND.ICHARM.GT.0) THEN
CZCZ    CALL ESDDERZ(ICHARM,X(LLAYFRG),X(LINDAT),X(LIABDFG),X(LJABDFG)
CZCZ *              ,X(LINDBD),X(LFMOZAN),X(LFMOC),ZNUC,CX,CY,CZ)
C       WRITE(6,*) 'WWWEXTZ',ICHARM,-ZNUC,CX,CY,CZ
CZCZ    IF(ABS(ZNUC).LT.1.0D-08) GOTO 2000
C       THE CHARGE IS OF COURSE EITHER 0 OR NOT, NOT MERELY CLOSE TO 0.
CZ      ENDIF
C
C     ----- I SHELL
C
      DO 1400 II = 1,NSHELL
C
C           GO PARALLEL!
C
        IF(NXT .AND. GOPARR) THEN
           MINE = MINE + 1
           IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
           IF(NEXT.NE.MINE) GO TO 1400
        END IF
C
        IAT = KATOM(II)
        XI = C(1,IAT)
        YI = C(2,IAT)
        ZI = C(3,IAT)
        I1 = KSTART(II)
        I2 = I1+KNG(II)-1
        LIT = KTYPE(II)
        MINI = KMIN(II)
        MAXI = KMAX(II)
        LOCI = KLOC(II)-MINI
        LITDER = LIT + 1
C
C     ----- J SHELL
C
        DO 1300 JJ = 1,NSHELL
C
C           GO PARALLEL!
C
          IF((.NOT.NXT) .AND. GOPARR) THEN
             IPCOUNT = IPCOUNT + 1
             IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 1300
          END IF
C
          JAT = KATOM(JJ)
          XJ = C(1,JAT)
          YJ = C(2,JAT)
          ZJ = C(3,JAT)
          J1 = KSTART(JJ)
          J2 = J1+KNG(JJ)-1
          LJT = KTYPE(JJ)
          MINJ = KMIN(JJ)
          MAXJ = KMAX(JJ)
          LOCJ = KLOC(JJ)-MINJ
CZCZ      LJTMOD = LJT+2
          NROOTS = (LIT+LJT-1)/2 + 1
          RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
C
C     ----- I PRIMITIVE
C
        DO 1100 IG = I1,I2
          AI = EX(IG)
          ARRI = AI*RR
          AXI = AI*XI
          AYI = AI*YI
          AZI = AI*ZI
          CSI=CS(IG)
          CPI=CP(IG)
          CDI=CD(IG)
          CFI=CF(IG)
          CGI=CG(IG)
C
C     ----- J PRIMITIVE
C
          DO 1000 JG = J1,J2
            AJ = EX(JG)
            AA = AI+AJ
            AA1 = ONE/AA
            DUM = AJ*ARRI*AA1
            IF(DUM .GT. TOL) GO TO 1000
            FAC = EXP(-DUM)
            CSJ = CS(JG)
            CPJ = CP(JG)
            CDJ = CD(JG)
            CFJ = CF(JG)
            CGJ = CG(JG)
            AX = (AXI+AJ*XJ)*AA1
            AY = (AYI+AJ*YJ)*AA1
            AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
            IJ = 0
            DUM1 = ZERO
            DUM2 = DUM1
            DO 600 I=MINI,MAXI
              IF(I.EQ.1) DUM1=CSI*FAC
              IF(I.EQ.2) DUM1=CPI*FAC
              IF(I.EQ.5) DUM1=CDI*FAC
              IF(I.EQ.8.AND.NORM) DUM1=DUM1*SQRT3
              IF(I.EQ.11) DUM1=CFI*FAC
              IF(I.EQ.14.AND.NORM) DUM1=DUM1*SQRT5
              IF(I.EQ.20.AND.NORM) DUM1=DUM1*SQRT3
              IF(I.EQ.21) DUM1=CGI*FAC
              IF(I.EQ.24.AND.NORM) DUM1=DUM1*SQRT7
              IF(I.EQ.30.AND.NORM) DUM1=DUM1*SQRT5/SQRT3
              IF(I.EQ.33.AND.NORM) DUM1=DUM1*SQRT3
C
              DO 600 J = MINJ,MAXJ
                IF(J.EQ.1) DUM2=DUM1*CSJ
                IF(J.EQ.2) DUM2=DUM1*CPJ
                IF(J.EQ.5) DUM2=DUM1*CDJ
                IF(J.EQ.8.AND.NORM) DUM2=DUM2*SQRT3
                IF(J.EQ.11) DUM2=DUM1*CFJ
                IF(J.EQ.14.AND.NORM) DUM2=DUM2*SQRT5
                IF(J.EQ.20.AND.NORM) DUM2=DUM2*SQRT3
                IF(J.EQ.21) DUM2=DUM1*CGJ
                IF(J.EQ.24.AND.NORM) DUM2=DUM2*SQRT7
                IF(J.EQ.30.AND.NORM) DUM2=DUM2*SQRT5/SQRT3
                IF(J.EQ.33.AND.NORM) DUM2=DUM2*SQRT3
C
                IJ=IJ+1
                IF(RMETHOD.EQ.RESC) THEN
                  DEN=ONE
                  IF(LOCI+I.EQ.LOCJ+J) DEN=TWO
                ELSE
                  NN = IA(MAX0(LOCI+I,LOCJ+J))+MIN0(LOCI+I,LOCJ+J)
                  DEN = DENAB(NN)
                  DEN = DEN + DEN
                  IF(RMETHOD.EQ.ANESC) THEN
                    IF(MODUS.EQ.1) DEN=DEN*FSC2
                    IF(MODUS.EQ.2) DEN=-DEN
                  ENDIF
                ENDIF
                DIJ(IJ)=DUM2*DEN
  600         CONTINUE
C
C     -----  KINETIC ENERGY
C
cz    IF(ICHARM.GT.0) GO TO 700
C
cz    T = SQRT(AA1)
cz    X0 = AX
cz    Y0 = AY
cz    Z0 = AZ
cz    DO 660 J = 1,LJTMOD
cz      NJ =J
cz      DO 660 I = 1,LITDER
cz        NI = I
cz        CALL VINT
cz        XS(I,J)=XINT*T
cz        YS(I,J)=YINT*T
cz        ZS(I,J)=ZINT*T
cz660 CONTINUE
cz    CALL DTXYZ(XT,YT,ZT,XS,YS,ZS,LITDER,LJT,AJ)
cz    CALL DERI(DXS,DYS,DZS,XS,YS,ZS,LIT,LJT,AI)
cz    CALL DERI(DXT,DYT,DZT,XT,YT,ZT,LIT,LJT,AI)
cz    IJ=0
cz    DO 680 I=MINI,MAXI
cz      IX=IJX(I)
cz      IY=IJY(I)
cz      IZ=IJZ(I)
cz      DO 670 J=MINJ,MAXJ
cz        JX=IJX(J)
cz        JY=IJY(J)
cz        JZ=IJZ(J)
cz        DUMX=DXT(IX,JX)* YS(IY,JY)* ZS(IZ,JZ)
cz   1        +DXS(IX,JX)* YT(IY,JY)* ZS(IZ,JZ)
cz   2        +DXS(IX,JX)* YS(IY,JY)* ZT(IZ,JZ)
cz        DUMY= XT(IX,JX)*DYS(IY,JY)* ZS(IZ,JZ)
cz   1        + XS(IX,JX)*DYT(IY,JY)* ZS(IZ,JZ)
cz   2        + XS(IX,JX)*DYS(IY,JY)* ZT(IZ,JZ)
cz        DUMZ= XT(IX,JX)* YS(IY,JY)*DZS(IZ,JZ)
cz   1        + XS(IX,JX)* YT(IY,JY)*DZS(IZ,JZ)
cz   2        + XS(IX,JX)* YS(IY,JY)*DZT(IZ,JZ)
cz        IJ=IJ+1
cz        IF(RMETHOD.EQ.RESC) THEN
cz           NN=IA(MAX0(LOCI+I,LOCJ+J))+MIN0(LOCI+I,LOCJ+J)
cz           DUMINT=DIJ(IJ)
cz           DTINT(NN,1,IAT)=DTINT(NN,1,IAT)+DUMX*DUMINT
cz           DTINT(NN,2,IAT)=DTINT(NN,2,IAT)+DUMY*DUMINT
cz           DTINT(NN,3,IAT)=DTINT(NN,3,IAT)+DUMZ*DUMINT
cz        ELSE
cz           DE(1,IAT)=DE(1,IAT)+ DUMX*DIJ(IJ)
cz           DE(2,IAT)=DE(2,IAT)+ DUMY*DIJ(IJ)
cz           DE(3,IAT)=DE(3,IAT)+ DUMZ*DIJ(IJ)
cz        ENDIF
cz670     CONTINUE
cz680   CONTINUE
cz      IF(RMETHOD.EQ.ANESC.AND.MODUS.NE.0) GOTO 1000
C
C     ..... NUCLEAR ATTRACTION
C
CZ700 CONTINUE
      AAX = AA*AX
      AAY = AA*AY
      AAZ = AA*AZ
CZ    IF(ICHARM.EQ.0) THEN
CZ       MAXIC=NAT
CZ       IF(ISEPS) MAXIC=NAT+NPS
CZ    ELSE
CZ       MAXIC=1
CZ    END IF
C
      DO 790 IC = 1,MAXIC
CZ
CZ       FMOATFRG2 is a clone of FMOATFRG located in HSANDT
CZ
         CALL FMOATFRG2(IC,X(LINDAT),X(LINDATG),X(LIAGLOB),
     *                  X(LIALOC),X(LIABDFG),X(LJABDFG),
     *                  X(LINDBD),X(LFMOZAN),X(LFMOC),
     *                  X(LPOPMAT),KATOM(II),KATOM(JJ),
     *                  KFG,ZNUC,CX,CY,CZ)
         IF (ABS(ZNUC).LT.1.0D-08) GOTO 790
CZ
CZ       IF(ICHARM.EQ.0) THEN
CZ          IF(IC.LE.NAT) THEN
CZ             ZNUC = -ZAN(IC)
CZ             CX = C(1,IC)
CZ             CY = C(2,IC)
CZ             CZ = C(3,IC)
CZ          ELSE
CZ             IAZ=IATSP(IC-NAT)
CZ             ZNUC = -COSZAN(IC-NAT)
CZ             CX = CORZAN(1,IC-NAT)
CZ             CY = CORZAN(2,IC-NAT)
CZ             CZ = CORZAN(3,IC-NAT)
CZ          END IF
CZ       ELSE
CZ         IF(.NOT.ESD) THEN
CZ          ZNUC = -QCHM(ICHARM)
CZ          CX = XCHM(ICHARM)
CZ          CY = YCHM(ICHARM)
CZ          CZ = ZCHM(ICHARM)
CZ          ENDIF
C           FOR ESD THE VARIABLES ARE SET ABOVE IN ESDDERZ.
CZ       END IF
         XX = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
         IF(NROOTS.LE.3) CALL RT123
         IF(NROOTS.EQ.4) CALL ROOT4
         IF(NROOTS.EQ.5) CALL ROOT5
         DO 740 K = 1,NROOTS
            UU = AA*U(K)
            WW = W(K)*ZNUC
            TT = ONE/(AA+UU)
            T = SQRT(TT)
            X0 = (AAX+UU*CX)*TT
            Y0 = (AAY+UU*CY)*TT
            Z0 = (AAZ+UU*CZ)*TT
            DO 730 J = 1,LJT
               NJ = J
               DO 720 I = 1,LITDER
                  NI = I
                  CALL VINT
                  XV(I,J,K) = XINT
                  YV(I,J,K) = YINT
                  ZV(I,J,K) = ZINT*WW
  720          CONTINUE
  730       CONTINUE
            CALL DERI(DXV(1,1,K),DYV(1,1,K),DZV(1,1,K),
     *                 XV(1,1,K), YV(1,1,K), ZV(1,1,K),LIT,LJT,AI)
  740    CONTINUE
         IJ=0
         DO 780 I=MINI,MAXI
           IX=IJX(I)
           IY=IJY(I)
           IZ=IJZ(I)
           DO 780 J=MINJ,MAXJ
             JX=IJX(J)
             JY=IJY(J)
             JZ=IJZ(J)
             DUMX=ZERO
             DUMY=ZERO
             DUMZ=ZERO
             DO 770 K=1,NROOTS
               DUMX=DUMX+DXV(IX,JX,K)* YV(IY,JY,K)* ZV(IZ,JZ,K)
               DUMY=DUMY+ XV(IX,JX,K)*DYV(IY,JY,K)* ZV(IZ,JZ,K)
               DUMZ=DUMZ+ XV(IX,JX,K)* YV(IY,JY,K)*DZV(IZ,JZ,K)
  770        CONTINUE
             IJ=IJ+1
CZ 
CZCZCZ ???   IF((IC.GT.NAT).AND.(IAT.EQ.IAZ)) GOTO 780
CZ
             DUMINT=DIJ(IJ)*AA1*PI212
             IF(RMETHOD.EQ.RESC) THEN
                NN=IA(MAX0(LOCI+I,LOCJ+J))+MIN0(LOCI+I,LOCJ+J)
                DVINT(NN,1,IAT)=DVINT(NN,1,IAT)+DUMX*DUMINT
                DVINT(NN,2,IAT)=DVINT(NN,2,IAT)+DUMY*DUMINT
                DVINT(NN,3,IAT)=DVINT(NN,3,IAT)+DUMZ*DUMINT
C               IE. NO CHARMS NOW WITH RESC
             ELSE
CZ              IF(ESD.AND.ICHARM.GT.0) THEN
CZ              SAVE TVDER GRADIENT TERMS TO IBODY ARRAYS
                DDDDX = DUMX*DUMINT
                DDDDY = DUMY*DUMINT
                DDDDZ = DUMZ*DUMINT
                IAZ=IXFTCH(X(LIAGLOB),IAT)
                if (idamm.ne.0) then
                  damm = idamm
                  IOFF=LFMODE         +(IAZ-1)*3
                  X(IOFF  )=X(IOFF  )+damm*DDDDX
                  X(IOFF+1)=X(IOFF+1)+damm*DDDDY
                  X(IOFF+2)=X(IOFF+2)+damm*DDDDZ
                end if
                if (idadd.ne.0) then
                  dadd = idadd
                  IOFF=LFMODE+3*NATFMO+(IAZ-1)*3
                  X(IOFF  )=X(IOFF  )+dadd*DDDDX
                  X(IOFF+1)=X(IOFF+1)+dadd*DDDDY
                  X(IOFF+2)=X(IOFF+2)+dadd*DDDDZ
                end if
                if (idatt.ne.0) then
                  datt = idatt
                  IOFF=LFMODE+6*NATFMO+(IAZ-1)*3
                  X(IOFF  )=X(IOFF  )+datt*DDDDX
                  X(IOFF+1)=X(IOFF+1)+datt*DDDDY
                  X(IOFF+2)=X(IOFF+2)+datt*DDDDZ
                end if
CZCZ               IAZ=IXFTCH(X(LIAGLOB),IAT)
CZCZ               IOFF=LFMODE+3*NATFMO+(IAZ-1)*3
CZCZ               X(IOFF  )=X(IOFF  )+DUMX*DUMINT
CZCZ               X(IOFF+1)=X(IOFF+1)+DUMY*DUMINT
CZCZ               X(IOFF+2)=X(IOFF+2)+DUMZ*DUMINT
C                  WRITE(6,*) 'WWWADD1',IAZ,IAT
CZ              ELSE
CZ                 I wonder if this is tvder gradient
CZ                 DE(1,IAT)=DE(1,IAT)+DUMX*DUMINT
CZ                 DE(2,IAT)=DE(2,IAT)+DUMY*DUMINT
CZ                 DE(3,IAT)=DE(3,IAT)+DUMZ*DUMINT
CZ              ENDIF
CZCZ            IF(ISEPS.AND.(IC.GT.NAT)) THEN
CZCZ               DE(1,IAZ)=DE(1,IAZ)-DUMX*(DIJ(IJ)*AA1*PI212)
CZCZ               DE(2,IAZ)=DE(2,IAZ)-DUMY*(DIJ(IJ)*AA1*PI212)
CZCZ               DE(3,IAZ)=DE(3,IAZ)-DUMZ*(DIJ(IJ)*AA1*PI212)
CZCZ            END IF
CZ              IF(ICHARM.GT.0) THEN
CZ                 IF(ESD) THEN
CZ              SAVE HELFEY GRADIENT TERMS TO IBODY ARRAYS
                if (idamm.ne.0) then
                  damm = idamm
                  IOFF=LFMODE         +(IC-1)*3
                  X(IOFF  )=X(IOFF  )-damm*DDDDX
                  X(IOFF+1)=X(IOFF+1)-damm*DDDDY
                  X(IOFF+2)=X(IOFF+2)-damm*DDDDZ
                end if
                if (idadd.ne.0) then
                  dadd = idadd
                  IOFF=LFMODE+3*NATFMO+(IC-1)*3
                  X(IOFF  )=X(IOFF  )-dadd*DDDDX
                  X(IOFF+1)=X(IOFF+1)-dadd*DDDDY
                  X(IOFF+2)=X(IOFF+2)-dadd*DDDDZ
                end if
                if (idatt.ne.0) then
                  datt = idatt
                  IOFF=LFMODE+6*NATFMO+(IC-1)*3
                  X(IOFF  )=X(IOFF  )-datt*DDDDX
                  X(IOFF+1)=X(IOFF+1)-datt*DDDDY
                  X(IOFF+2)=X(IOFF+2)-datt*DDDDZ
                end if
CZCZ                  IOFF=LFMODE+3*NATFMO+(ICHARM-1)*3
CZCZ                  X(IOFF  )=X(IOFF  )-DUMX*DUMINT
CZCZ                  X(IOFF+1)=X(IOFF+1)-DUMY*DUMINT
CZCZ                  X(IOFF+2)=X(IOFF+2)-DUMZ*DUMINT
C                     WRITE(6,*) 'WWWADD2',ICHARM
C            CHARMM CODE BELOW NEEDS THE FORCE, NOT THE GRADIENT?
CZ                 ELSE
CZ                    DXELMM(ICHARM)=DXELMM(ICHARM)+DUMX*DUMINT
CZ                    DYELMM(ICHARM)=DYELMM(ICHARM)+DUMY*DUMINT
CZ                    DZELMM(ICHARM)=DZELMM(ICHARM)+DUMZ*DUMINT
CZ                 ENDIF
CZ             END IF
             ENDIF
  780    CONTINUE
  790 CONTINUE
C
 1000 CONTINUE
 1100 CONTINUE
C
C     ----- END OF PRIMITIVE LOOPS -----
C
      IF(DBG) THEN
         WRITE(IW,9100) II,JJ
         CALL EGOUT(DE,NAT)
      END IF
 1300 CONTINUE
 1400 CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
CZCZ
C2000 CONTINUE
CZCZ
      IF(OUT) THEN
         IF(GOPARR) CALL DDI_GSUMF(1503,DE,3*NAT)
         IF(MASWRK) THEN
            WRITE(IW,9000)
            CALL EGOUT(DE,NAT)
         END IF
         IF(GOPARR) CALL DSCAL(3*NAT,ONE/NPROC,DE,1)
C              WILL ONLY PRINT PARTIAL CONTRIBUTIONS IN PARALLEL.
CZCZ     IF(ESD) THEN
         IF(.true.) THEN
            IF(MASWRK) WRITE(IW,9088)
            CALL EGOUT(X(LFMODE+3*NATFMO),NATFMO)
         END IF
      END IF
CZCZ     IF(ESD) THEN
CZCZ for debug
         DEBUGOPT=.FALSE.
         IF(DEBUGOPT) THEN
            IF(MASWRK) WRITE(IW,9088)
            IF(MASWRK) WRITE(IW,*) 'NATNAT',NAT
            CALL EGOUT(X(LFMODE),NATFMO*nbody)
         END IF
CZCZ
C
CZ
      IF(GOPARR .AND. NXT) CALL DDI_DLBRESET
C
      IF(RMETHOD.EQ.RESC) THEN
        IF(GOPARR) CALL DDI_GSUMF(1503,DTINT,L2*3*NAT)
        IF(GOPARR) CALL DDI_GSUMF(1504,DVINT,L2*3*NAT)
        SVDSKW=DSKWRK
        DSKWRK=.TRUE.
         LREC=3*NAT
C        STARTING RECORD FOR THE T DERIVATIVE INTEGRALS
        DO I=1,NAT
          DO J=1,3
            LREC=LREC+1
            CALL RAWRIT(NQRDAF,X(MORDA),DTINT(1,J,I),L2,LREC,0)
C
C           ADD TWO CONTRIBUTIONS, HELLMANN-FEYNMAN AND DV GOT HERE
C
            CALL RAREAD(NQRDAF,X(MORDA),DTINT(1,J,I),L2,LREC+6*NAT,0)
            CALL DAXPY(L2,ONE,DTINT(1,J,I),1,DVINT(1,J,I),1)
            CALL RAWRIT(NQRDAF,X(MORDA),DVINT(1,J,I),L2,LREC+6*NAT,0)
          ENDDO
        ENDDO
        DSKWRK=SVDSKW
      ENDIF
C
      RETURN
C
 9000 FORMAT(/10X,38(1H-)/10X,'GRADIENT INCLUDING AO DERIVATIVE TERMS'/
     *        10X,38(1H-))
 9088 FORMAT(/10X,36(1H-)/10X,'ONE-ELECTRON ESD GRADIENT',/10X,36(1H-))
 9100 FORMAT(1X,'TVDER: SHELLS II,JJ=',2I5)
      END
C*MODULE FMOESD  *DECK FMOATFRG2
CZ    SUBROUTINE fmoatfrg(iat0,indat,indatg,iaglob,ialoc,iabdfg,jabdfg,
CZ   *                    indbd,fmozan,fmoc,popmat,iats,jats,fracv,
CZ   *                    fracesp,iatfrg,ZNUC,cx,cy,cz)
      SUBROUTINE fmoatfrg2(iat0,indat,indatg,iaglob,ialoc,iabdfg,jabdfg,
     *                    indbd,fmozan,fmoc,popmat,iats,jats,iatfrg,
     *                    ZNUC,cx,cy,cz)
CZ
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      logical iainside,jainside,smartr(2),bimer(3)
      dimension indat(*),indatg(natfmo,*),iaglob(*),iabdfg(*),jabdfg(*),
     *          indbd(maxabd,*),fmozan(*),fmoc(3,*),popmat(maxnat,nfg,2)
     *         ,ialoc(*)
      Parameter (zero=0.0D+00,one=1.0D+00)
C
      COMMON /FMCOM / XX(1)
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
      COMMON /MCPFMO/ IMCPFMO,LFZCOR,LIFMPTYP,LIFMPTYP2,LMCPSW,
     *                LIZCOR2
C
      data bimer/3*.false./
CZ
CZ    This is a clone of fmoatfrg
CZ
c
c     Determine if atom IAT belongs to fragment 
c
      if(ifmostp.eq.2.and.iand(ixesp,4096).ne.0) then
        iatfrg=indat(iat0)
        rk=fmodist(icurfg,0,0,iatfrg)
        if(icurfg.ne.iatfrg) then
CZCZ      fracv=0
          fracesp=0
          igot=0
          if(rk.eq.0) fracesp=one
          znuc=-fmozan(iat0)
          cx=fmoc(1,iat0)
          cy=fmoc(2,iat0)
          cz=fmoc(3,iat0)
          if(nbdfg.eq.0) goto 90
          do iabd=1,maxabd
            ibdfg=indbd(iabd,iat0)
            if(ibdfg.eq.0) goto 90
            ia=abs(iabdfg(ibdfg))
            ja=abs(jabdfg(ibdfg))
            iafrg=indat(ia)
            jafrg=indat(ja)
            if(iafrg.eq.icurfg.or.jafrg.eq.icurfg) then
c             this is handled below
              igot=1
              goto 90
            endif
c           ri=fmodist(icurfg,0,0,iafrg)
            rj=fmodist(icurfg,0,0,jafrg)
c           The code below may not work for multiply cut BDAs?
            if(rj.eq.0) then
              if(rk.ne.0) fracesp=fracesp-one/znuc
            else
              if(rk.eq.0) fracesp=fracesp+one/znuc
            endif
          enddo
   90     continue
          if(igot.eq.0) return
        endif 
      endif 
CZ    if(ifmostp.ne.6) then
        iat=iat0
CZCZ    znuc=-fmozan(iat)
        IF (IMCPFMO.EQ.1) THEN
          znuc=-(fmozan(iat)-XX(LFZCOR+iat-1))
c         write(6,*) 'ZNUC', znuc, fmozan(iat)
        ELSE
          znuc=-fmozan(iat)
        END IF
CZ
        cx=fmoc(1,iat)
        cy=fmoc(2,iat)
        cz=fmoc(3,iat)
        iz=1
        ifg=icurfg
        factk=one
CZ    else
CZ      iat=iat0+nat
CZ      znuc=-zan(iat)
CZ      cx=c(1,iat)
CZ      cy=c(2,iat)
CZ      cz=c(3,iat)
c       iat below must be global iat
CZ      iat=iaglob(iat)
CZ      iz=2
CZ      ifg=ncursh
CZ      factk=half
CZ    endif  
      jfg=jcurfg
      lfg=kcurfg
      iatfrg=indat(iat)
      smartr(1)=iand(modesp,7).eq.1.and.jfg.ne.0
      smartr(2)=iand(modesp,7).eq.2.and.jfg.ne.0
      if(smartr(1).and.nbdfg.ne.0) then
        bimer(1)=fmodist(ifg,0,0,jfg).eq.0
        if(lfg.eq.0) then
          if(bimer(1)) smartr(1)=.false.
        else
          bimer(2)=fmodist(ifg,0,0,lfg).eq.0
          bimer(3)=fmodist(jfg,0,0,lfg).eq.0
        if(bimer(1).and.(bimer(2).or.bimer(3)).or.bimer(2).and.bimer(3))
     *    smartr(1)=.false.
        endif 
      endif 
      iatsg=iaglob(iats)
      jatsg=iaglob(jats)
      iifg=indat(iatsg)
      jjfg=indat(jatsg)
c     3 -body terms do not work with BSSE.
CZCZ  if(nbdfg.eq.0.or.ifmostp.eq.6) then
      if(nbdfg.eq.0) then
CZ      fracv=zero
        fracesp=one
c       BSSE
        if(ifmostp.eq.5) then
          if(nbsse.eq.1.and.iatfrg.eq.jfg .or. nbsse.eq.2.and.
     *      iatfrg.ne.ifg.and.(iifg.eq.jfg.or.jjfg.eq.jfg)) fracesp=zero
        endif
c       decide if the point charge iat should be included into ESP (keep=1)
CZ      if(ifmostp.ne.6) then
         if(fracesp.eq.one.and.(nbsse.ne.2.or.ifmostp.ne.5.or.iatfrg.ne.
     *  jfg).and.(iatfrg.eq.ifg.or.iatfrg.eq.jfg.or.iatfrg.eq.lfg)) then
CZ        fracv=one
          fracesp=zero
         endif 
CZ      endif
      else
c       this code does not support nbsse=1,2 and 3.
        if(nbsse.ne.0) then
          write(6,*) 'bad nbsse',nbsse
          call abrt
        endif
c       set values for the case atom IAT is not shared
        if(iatfrg.eq.ifg.or.iatfrg.eq.jfg.or.iatfrg.eq.lfg) then
CZ        fracv=one
          fracesp=zero
        else
CZ        fracv=zero
          fracesp=one
        endif
        fracesp1=zero
c       do ibdfg=1,nbdfg
c       loop over all broken bonds in which atom IAT is involved.
c       we are only interested really in one n-mer here (n=1,2), so
c       we find if an atom is split between this n-mer and the rest, then quit. 
        do iabd=1,maxabd
          ibdfg=indbd(iabd,iat)
          if(ibdfg.eq.0) goto 100
          ia=abs(iabdfg(ibdfg))
          ja=abs(jabdfg(ibdfg))
c         reordering is now done in fmobon
c         the code below assumes the canonical order (left is negative). 
c         if(ja.lt.0) then
c           if(ia.lt.0) call abrt
c           idum=ia
c           ia=ja
c           ja=idum
c         endif
c         ia=abs(ia)
          iafrg=indat(ia)
          jafrg=indat(ja)
          iainside=iafrg.eq.ifg.or.iafrg.eq.jfg.or.iafrg.eq.lfg
          jainside=jafrg.eq.ifg.or.jafrg.eq.jfg.or.jafrg.eq.lfg
c         exclude the case when a cut bond is inside a dimer
          if(iat.eq.ia.and..not.(iainside.and.jainside)) then
            if(iainside) then
CZ            fracv=(znuc+one)/znuc
              fracesp=-one/znuc
            endif
            if(jainside) then
CZ            fracv=-one/znuc
              fracesp=(znuc+one)/znuc
            endif
          endif
c
c           now add partial charges. 
c           Here we add the "right" side (that is, pseudoproton). 
c           The left side (N-1) is treated below, along with all other cases. 
c           
          if(iat.eq.ia.and..not.jainside.and.resppc(iz).ne.zero) then
c           kfg is equal to jafrg.
c           tricky part! kat comes from the ghost atom stored in ialoc
c           Note that this atom has coordinates equal to those of iat
c           but its charge is stored as processed below.
            kfg=ishft(ialoc(natfmo+ibdfg),-16)
            kat=iand(ialoc(natfmo+ibdfg),65535)
cnb         ifmostp.eq.6 cannot come here?!
            if(ifmostp.ne.6) then
              if(smartr(1).or.smartr(2)) then
                rk=fmosdist(iifg,jjfg,indatg(iatsg,1),indatg(jatsg,1),
     *                      ifg,jfg,lfg,kfg,bimer)
              else
                rk=fmodist(ifg,jfg,lfg,kfg)
              endif
            else
              rk=fmodist(ifg,0,0,kfg)
            endif
            if(rk.gt.resppc(iz)) then 
              fracesp1=fracesp1+popmat(kat,kfg,icurpop)*factk/znuc
            endif 
c           goto 100
          endif
          if(ifmostp.eq.6.and.resppc(iz).ne.zero) then
            write(6,*) 'fmoptc is not programmed for resppc(2) yet'
            call abrt
          endif
        enddo
  100   continue
        fracesp=fracesp+fracesp1
      endif
c     Add approximate 2e ESP contributions that become atomic Mulliken charges
c     centred at the atomic coordinates (point charge approximation).
      if(resppc(iz).ne.zero) then
        ind=iat
        kfg=ishft(ialoc(ind),-16)
        kat=iand(ialoc(ind),65535)
c       check if the gun is accidently or mistakenly loaded.  
        if(kfg.eq.0) call abrt 
c         
c       for regular runs we want all charges except those from I and J. 
CZCZ    if(ifmostp.ne.6.and.(ifg.eq.kfg.or.jfg.eq.kfg.or.lfg.eq.kfg)) 
CZCZ *  goto 200
        if(ifg.eq.kfg.or.jfg.eq.kfg.or.lfg.eq.kfg) goto 200
CZ
c       for esdim we only want charges coming from J
c
        if(ifmostp.ne.6) then
          if(smartr(1).or.smartr(2)) then
            rk=fmosdist(iifg,jjfg,indatg(iatsg,1),indatg(jatsg,1),
     *                   ifg,jfg,lfg,kfg,bimer)
          else
            rk=fmodist(ifg,jfg,lfg,kfg)
          endif
        else
          rk=fmodist(ifg,0,0,kfg)
        endif
        if(rk.gt.resppc(iz)) then 
          fracesp=fracesp+popmat(kat,kfg,icurpop)*factk/znuc
        endif 
  200   continue
      endif
CZ
CZ    FRACESP IS NOT NEEDED FOR ES1DER
CZ
      znuc = znuc * fracesp
CZ
      RETURN
      END
C*MODULE FMOESD  *DECK EXTSDER
      SUBROUTINE EXTSDER(EPS,DSINT,L1,L2,NOCON,LMONLAG,num0)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DBG,NORM,OUT
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT,SVDSKW,NOCON
CZ
      LOGICAL LMONLAG
CZ
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192)
C
      DIMENSION EPS(L2),DSINT(L2,3,NAT)
      DIMENSION DIJ(225), DIJINT(225), IJX(35), IJY(35),IJZ(35),
     *           XS(6,5), YS(6,5), ZS(6,5),
     *          DXS(5,5),DYS(5,5),DZS(5,5)
C
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /FMCOM / X(1)
CZ
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
CZ
C
      PARAMETER (SQRT3=1.73205080756888D+00)
      PARAMETER (SQRT5=2.23606797749979D+00, SQRT7=2.64575131106459D+00)
      PARAMETER (ONE=1.0D+00)
      PARAMETER (RLN10=2.30258D+00)
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HSDER    /, GRD1/8HGRD1    /,
     *     RESC/8HRESC    /
C
      DATA IJX / 1, 2, 1, 1, 3, 1, 1, 2, 2, 1,
     *           4, 1, 1, 3, 3, 2, 1, 2, 1, 2,
     *           5, 1, 1, 4, 4, 2, 1, 2, 1, 3,
     *           3, 1, 3, 2, 2/
      DATA IJY / 1, 1, 2, 1, 1, 3, 1, 2, 1, 2,
     *           1, 4, 1, 2, 1, 3, 3, 1, 2, 2,
     *           1, 5, 1, 2, 1, 4, 4, 1, 2, 3,
     *           1, 3, 2, 3, 2/
      DATA IJZ / 1, 1, 1, 2, 1, 1, 3, 1, 2, 2,
     *           1, 1, 4, 1, 2, 1, 2, 3, 3, 2,
     *           1, 1, 5, 1, 2, 1, 2, 4, 4, 1,
     *           3, 3, 2, 2, 3/
C
C     ----- CALCULATE DERIVATIVES OF THE OVERLAP MATRIX -----
C     THIS TERM RESULTS FROM DERIVATIVES OF THE DENSITY MATRIX,
C     WHICH CAN IN THE CASE OF A SYMMETRIC LAGRANGIAN MATRIX, CAN
C     BE CONVERTED INSTEAD INTO DERIVATIVES OF OVERLAP INTEGRALS.
C     RESC RUNS NEED BOTH CONTRACTION OF THE INTEGRALS WITH DENSITY
C     AND THE INTEGRALS THEMSELVES.
C     IN CASE OF INTERNALLY UNCONTRACTED RESC, ONE HAS TO TRANSFORM
C     THE OVERLAP DERIVATIVES TO CONTRACTED BASIS BEFORE
C     SUMMING (CONTRACTING) OVERLAP DERIVATIVES WITH THE LAGRANGIAN
C
C     NOCON IS TRUE IF LAGRANGIAN IS NOT TO BE SUMMED WITH.
C
      DBG = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME)  .AND.  MASWRK
      OUT = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME  .OR.
     *       EXETYP.EQ.GRD1   .OR.  NPRINT.EQ.-3)
      IF(DBG) WRITE(IW,9000)
C
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
CZCZ  IF(RMETHOD.EQ.RESC) CALL VCLR(DSINT,1,L2*3*NAT)
      IF(LMONLAG) CALL VCLR(DSINT,1,L2*3*NAT)
CZ
C
C     INITIALIZE PARALLEL
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT = -1
      MINE = -1
C
      IF(DBG) THEN
         WRITE(IW,9010)
         CALL PRTRI(EPS,L1)
      END IF
C
C     ----- I SHELL
C
CZCZ  DO 780 II = 1,NSHELL
      DO 780 II = NCURSH+1,NSHELL
CZ
C
C           GO PARALLEL!
C
      IF(NXT .AND. GOPARR) THEN
         MINE = MINE + 1
         IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
         IF(NEXT.NE.MINE) GO TO 780
      END IF
C
      IAT = KATOM(II)
      XI = C(1,IAT)
      YI = C(2,IAT)
      ZI = C(3,IAT)
      I1 = KSTART(II)
      I2 = I1+KNG(II)-1
      LIT = KTYPE(II)
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II)-MINI
      LITDER = LIT+1
C
C     ----- J SHELL
C
CZCZ  DO 760 JJ = 1,II
      DO 760 JJ = NCURSH+1,II
CZ
C
C           GO PARALLEL!
C
      IF((.NOT.NXT) .AND. GOPARR) THEN
         IPCOUNT = IPCOUNT + 1
         IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 760
      END IF
C
      JAT = KATOM(JJ)
      XJ = C(1,JAT)
      YJ = C(2,JAT)
      ZJ = C(3,JAT)
      J1 = KSTART(JJ)
      J2 = J1+KNG(JJ)-1
      LJT = KTYPE(JJ)
      MINJ = KMIN(JJ)
      MAXJ = KMAX(JJ)
      LOCJ = KLOC(JJ)-MINJ
      RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
      IF(II.EQ.JJ) GO TO 760
C
C     ----- I PRIMITIVE
C
      DO 640 IG = I1,I2
        AI = EX(IG)
        ARRI = AI*RR
        AXI = AI*XI
        AYI = AI*YI
        AZI = AI*ZI
        CSI=CS(IG)
        CPI=CP(IG)
        CDI=CD(IG)
        CFI=CF(IG)
        CGI=CG(IG)
C
C     ----- J PRIMITIVE
C
        DO 620 JG = J1,J2
          AJ = EX(JG)
          AA = AI+AJ
          AA1 = ONE/AA
          DUM = AJ*ARRI*AA1
          IF(DUM .GT. TOL) GO TO 620
          FAC = EXP(-DUM)
          CSJ = CS(JG)
          CPJ = CP(JG)
          CDJ = CD(JG)
          CFJ = CF(JG)
          CGJ = CG(JG)
          AX = (AXI+AJ*XJ)*AA1
          AY = (AYI+AJ*YJ)*AA1
          AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
          IJ = 0
          DO 580 I = MINI,MAXI
            IF(I.EQ.1) DUM1=CSI*FAC
            IF(I.EQ.2) DUM1=CPI*FAC
            IF(I.EQ.5) DUM1=CDI*FAC
            IF(I.EQ.8.AND.NORM) DUM1=DUM1*SQRT3
            IF(I.EQ.11) DUM1=CFI*FAC
            IF(I.EQ.14.AND.NORM) DUM1=DUM1*SQRT5
            IF(I.EQ.20.AND.NORM) DUM1=DUM1*SQRT3
            IF(I.EQ.21) DUM1=CGI*FAC
            IF(I.EQ.24.AND.NORM) DUM1=DUM1*SQRT7
            IF(I.EQ.30.AND.NORM) DUM1=DUM1*SQRT5/SQRT3
            IF(I.EQ.33.AND.NORM) DUM1=DUM1*SQRT3
            DO 570 J = MINJ,MAXJ
              IF(J.EQ.1) DUM2=DUM1*CSJ
              IF(J.EQ.2) DUM2=DUM1*CPJ
              IF(J.EQ.5) DUM2=DUM1*CDJ
              IF(J.EQ.8.AND.NORM) DUM2=DUM2*SQRT3
              IF(J.EQ.11) DUM2=DUM1*CFJ
              IF(J.EQ.14.AND.NORM) DUM2=DUM2*SQRT5
              IF(J.EQ.20.AND.NORM) DUM2=DUM2*SQRT3
              IF(J.EQ.21) DUM2=DUM1*CGJ
              IF(J.EQ.24.AND.NORM) DUM2=DUM2*SQRT7
              IF(J.EQ.30.AND.NORM) DUM2=DUM2*SQRT5/SQRT3
              IF(J.EQ.33.AND.NORM) DUM2=DUM2*SQRT3
              IJ=IJ+1
CZCZ          NN=IA(LOCI+I)+(LOCJ+J)
              NN=IA(LOCI+I-NUM0)+(LOCJ+J-NUM0)
CZCZ djust the position of EPS()
CZCZ          IF(RMETHOD.EQ.RESC) THEN
              IF(LMONLAG) THEN
                DIJINT(IJ)=DUM2
              ENDIF
              DEN = EPS(NN)
              DEN=DEN+DEN
              DIJ(IJ) = DUM2*DEN
  570       CONTINUE
  580     CONTINUE
C
C     ----- OVERLAP
C
      T = SQRT(AA1)
      X0 = AX
      Y0 = AY
      Z0 = AZ
      DO 590 J = 1,LJT
        NJ = J
        DO 590 I = 1,LITDER
          NI = I
          CALL VINT
          XS(I,J)=XINT*T
          YS(I,J)=YINT*T
          ZS(I,J)=ZINT*T
  590 CONTINUE
C
      CALL DERI(DXS,DYS,DZS,XS,YS,ZS,LIT,LJT,AI)
C
      IJ=0
      DO 600 I=MINI,MAXI
        IX=IJX(I)
        IY=IJY(I)
        IZ=IJZ(I)
        DO 600 J=MINJ,MAXJ
          JX=IJX(J)
          JY=IJY(J)
          JZ=IJZ(J)
          DUMX=DXS(IX,JX)* YS(IY,JY)* ZS(IZ,JZ)
          DUMY= XS(IX,JX)*DYS(IY,JY)* ZS(IZ,JZ)
          DUMZ= XS(IX,JX)* YS(IY,JY)*DZS(IZ,JZ)
          IJ=IJ+1
CZCZ      IF(RMETHOD.EQ.RESC) THEN
          IF(LMONLAG) THEN
            NN=IA(LOCI+I)+(LOCJ+J)
            DUMINT=DIJINT(IJ)
            DSINT(NN,1,IAT)=DSINT(NN,1,IAT)+DUMX*DUMINT
            DSINT(NN,2,IAT)=DSINT(NN,2,IAT)+DUMY*DUMINT
            DSINT(NN,3,IAT)=DSINT(NN,3,IAT)+DUMZ*DUMINT
            DSINT(NN,1,JAT)=DSINT(NN,1,JAT)-DUMX*DUMINT
            DSINT(NN,2,JAT)=DSINT(NN,2,JAT)-DUMY*DUMINT
            DSINT(NN,3,JAT)=DSINT(NN,3,JAT)-DUMZ*DUMINT
          ENDIF
          IF(.NOT.NOCON) THEN
            DE(1,IAT)=DE(1,IAT)+(DUMX*DIJ(IJ))
            DE(2,IAT)=DE(2,IAT)+(DUMY*DIJ(IJ))
            DE(3,IAT)=DE(3,IAT)+(DUMZ*DIJ(IJ))
            DE(1,JAT)=DE(1,JAT)-(DUMX*DIJ(IJ))
            DE(2,JAT)=DE(2,JAT)-(DUMY*DIJ(IJ))
            DE(3,JAT)=DE(3,JAT)-(DUMZ*DIJ(IJ))
          ENDIF
  600 CONTINUE
C
  620 CONTINUE
  640 CONTINUE
C
C     ----- END OF PRIMITIVE LOOPS -----
C
      IF(DBG) THEN
         WRITE(IW,9100) II,JJ
         CALL EGOUT(DE,NAT)
      END IF
  760 CONTINUE
  780 CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
      IF(OUT) THEN
         IF(GOPARR) CALL DDI_GSUMF(1501,DE,3*NAT)
         IF(MASWRK) THEN
            WRITE(IW,9000)
            CALL EGOUT(DE,NAT)
         END IF
         IF(GOPARR) CALL DSCAL(3*NAT,ONE/NPROC,DE,1)
      END IF
      IF(GOPARR .AND. NXT) CALL DDI_DLBRESET
C
      IF(RMETHOD.EQ.RESC) THEN
        IF(GOPARR) CALL DDI_GSUMF(1501,DSINT,L2*3*NAT)
        SVDSKW=DSKWRK
        DSKWRK=.TRUE.
        LREC=0
C       STARTING RECORD FOR THE S DERIVATIVE INTEGRALS
        DO I=1,NAT
          DO J=1,3
            LREC=LREC+1
            CALL RAWRIT(NQRDAF,X(MORDA),DSINT(1,J,I),L2,LREC,0)
          ENDDO
        ENDDO
        DSKWRK=SVDSKW
      ENDIF
C
      RETURN
C
 9000 FORMAT(/10X,33(1H-)/10X,'GRADIENT INCLUDING DENSITY FORCES'/
     *        10X,33(1H-))
 9010 FORMAT(1X,'THE LAGRANGIAN IN THE AO BASIS IS')
 9100 FORMAT(1X,'SDER: SHELLS II,JJ=',2I5)
      END
C*MODULE FMOESD  *DECK ESDNSDER
      SUBROUTINE ESDNSDER(l1,l2,da,db,dc,wrk,wrk1,layfrg,idmrec0,
     *                    ESPA,WRKESP,numfrg,
     *                    mapi,mapj,map3,maxl1d,irec0,
     *                    orbxch1,enexch1,kodexch,jodexch,iodexch)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      Parameter (MXATM=2000,one=1.0D+00,half=0.5D+00)
C
      logical GOPARR,DSKWRK,MASWRK,doddcor
      logical orbxch1,enexch1,kodexch,jodexch,iodexch,some1
C
      dimension layfrg(*),idmrec0(*)
C
      DIMENSION DA(1),DB(1),DC(1)
      DIMENSION ESPA(1),WRK(1),WRK1(1),WRKESP(1)
      DIMENSION NUMFRG(1)
      DIMENSION mapi(1),mapj(1),map3(maxl1d,3) 
C
      COMMON /FMCOM / X(1)
      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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      common /fmopnt/ lichfg,lmulfg,lidmrec,lfrgnam,llayfrg,lindat,
     *                lnCBS,lfmozan,lfmoc,lfmomas,lizbas,liaglob,libdgh,
     *                liabdfg,ljabdfg,lnCAO,lidxCAO,liaprjo,ljaprjo,
     *                lCoreAO,lOccCor,lshiftb,liodfmo,lfmoda,lfmodb,
     *                lfmoespa,lfmoespb,llocfmo,lscffrg,lrij,lpopmul,
     *                lpopmat,lialoc,lindbd,liatfrg,lindfrg,lindgfrg,
     *                lnatfrg,lnat0frg,lianfrg,lzanfrg,lcfrg,llibish,
     *                llibnsh,llibng,lindatg,lfmobuf(3),lfmode,lnumfrg,
     *                lloctat,liaoglob,lloadm,lfmoge,ldgrid,liodcfmo,
     *                ljob2grp,lfmopg,lemocdr,luntxyz,luntrot,lstonep,
     *                lmapsu,lfrgmul
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
CZ
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
CZ
CZ    DENSITY DERIVATIVES OF EXTERNAL ELECTROSTATIC POTENTIALS
CZ
      some1 = .false.
      ifg  = icurfg
      jfg  = jcurfg
      kfg  = kcurfg
      ilay = icurlay
      L1I  = IAND(NUMFRG(IFG),65535)
      NAI  = ISHFT(NUMFRG(IFG),-16)
      NATI = ixftch(x(lnatfrg),ifg)
      IF (jfg.ne.0) THEN
        L1J  = IAND(NUMFRG(JFG),65535)
        NAJ  = ISHFT(NUMFRG(JFG),-16)
        NATJ = ixftch(x(lnatfrg),jfg)
      END IF
      IF (kfg.ne.0) THEN
        L1K  = IAND(NUMFRG(KFG),65535)
        NAK  = ISHFT(NUMFRG(KFG),-16)
        NATK = ixftch(x(lnatfrg),kfg)
      END IF

      IF (MPLEVL.EQ.0) THEN
        CALL DAREAD(IDAF,IODA,WRKESP,L2,16,0)
      ELSE
        CALL DAREAD(IDAF,IODA,WRKESP,L2,308,0)
      END IF
CZ
CZ    CONSTRUCT DI + DJ (+ DK) BLOCK DIAGONAL MATRIX (STORED IN ESPB) 
CZ
      IF (IFMOSTP.EQ.4.OR.IFMOSTP.EQ.8) THEN
        doddcor=iand(ixesp,1024).ne.0
        irec0 = idmrec0(ilay)
        if(doddcor) then
          irec0rhf=1
          if(irec0.eq.1) irec0rhf=nfg*2+1
          irec0=irec0rhf
        endif
      END IF

      NN = 1 
      IF (IFMOSTP.EQ.4) THEN
        NN = 2
        call diminid(0,jfg,ifg,l1j,l1i,naj,nai,l1,DA,db,wrk,wrk1,
     *               x(liodfmo),irec0,x(liabdfg),x(ljabdfg),x(lindat),
     *               x(liaglob),x(llocfmo),mapj,mapi,
     *               orbxch1,enexch1,jodexch,iodexch,0,dum,some1)
      ELSEIF (IFMOSTP.EQ.8) THEN
        NN = 3
        call triminid(kfg,jfg,ifg,l1k,l1j,l1i,nak,naj,nai,l1,DA,db,dc,
     *                wrk,wrk1,x(liodfmo),irec0,x(liabdfg),x(ljabdfg),
     *                x(lindat),x(liaglob),x(llocfmo),
     *                map3(1,3),map3(1,2),map3(1,1),
     *                orbxch1,enexch1,kodexch,jodexch,iodexch,some1)
      END IF
CZ
CZ    note that diminid and triminid destroy NO. 16 of Dictionary file
CZ

CZ
CZ    CALCULATE -(DI*VIJ*DI+DJ*VIJ*DJ)*SDER FOR DIMER
CZ
      CALL VCLR(DE,1,NAT*3)
CZ
CZ    CONSTRUCT MONOMER LAGRANGIAN TERM
CZ    NOTE: since the terms relative to fmoespa are calculated here,
CZ    the terms must be subtracted using FMOLAG in grad1.src. 
CZ
c     IF (NBDFG.NE.0) THEN
        IF (IFMOSTP.EQ.4) THEN
          natf2 = 3*natfmo
          CALL ESLAGRN(l1,l1i,l1j,l1j,nai,naj,nak,NATI,NATJ,NATJ,
     *                 orbxch1,iodexch,jodexch,kodexch,wrk,db,
     *                 layfrg,idmrec0,ESPA,
     *                 mapi,mapj,mapi,
     *                 irec0,x(liaglob),x(lfmode+natf2))
        ELSEIF (IFMOSTP.EQ.8) THEN
          natf3 = 6*natfmo
          CALL ESLAGRN(l1,l1i,l1j,l1k,nai,naj,nak,NATI,NATJ,NATK,
     *                 orbxch1,iodexch,jodexch,kodexch,wrk,db,
     *                 layfrg,idmrec0,ESPA,
     *                 map3(1,1),map3(1,2),map3(1,3),
     *                 irec0,x(liaglob),x(lfmode+natf3))
        END IF
c     ELSE
c       IF (NN.GE.2) THEN
c         CALL CPYTSQ(DA,DA(L2+1),L1,1)
c         CALL TFTRI(WRK,ESPA,DA(L2+1),WRK(L2+1),L1,L1,L1)
c         CALL DSCAL(L2, HALF,WRK,1)
c         l2i = (l1i*l1i+l1i)/2
c         do i = 1, L1j
c           kk = l2i+l1i*(i-1) + i*(i-1)/2
c           do j = 1, L1i
c             wrk(kk+j) = zero
c           end do 
c         end do
c         CALL SDER(WRK,WRK(L2+1),L1,L2,.FALSE.)
c       END IF
c     END IF
CZ
CZ    DI, DIJ OR DIJK DERIVATIVES 
CZ
c     if (maswrk) then
c       write (6,*) 'density2',L1,L1I,L1J
c       do i = 1, L1
c         k = (i*i-i)/2
c         write (6,'(6F12.8)') (da(k+j),j=1,i) 
c       end do
c     end if
      CALL CPYTSQ(WRKESP,DA(L2+1),L1,1)
      CALL TFTRI(WRK,ESPA,DA(L2+1),WRK(L2+1),L1,L1,L1)
      CALL DSCAL(L2,-HALF,WRK,1)
      CALL SDER(WRK,WRK(L2+1),L1,L2,.FALSE.)

      IF(GOPARR) CALL DDI_GSUMF(2418,DE,NAT*3)

      IF (MASWRK) THEN
        IF (NN.EQ.1) THEN
          CALL FMODEG(1,x(LFMODE),x(lfmopg),x(liaglob))
        ELSEIF (NN.EQ.2) THEN
          NFMOPT = LFMODE + NATFMO*3
          CALL FMODEG(1,x(NFMOPT),x(lfmopg),x(liaglob))
          NFMOPT = NFMOPT + NATFMO*3
          IF (NBODY.GT.2) CALL FMODEG(0,x(NFMOPT),x(lfmopg),x(liaglob))
          ! above must be changed
        ELSEIF (NN.EQ.3) THEN
          NFMOPT = LFMODE + NATFMO*6
          CALL FMODEG(1,x(NFMOPT),x(lfmopg),x(liaglob))
        END IF
      END IF

CZ
CZ    NOW DA IS DELTA DIJ WHICH IS USED FOR OTHER DERIVATIVES
CZ
      IF (NN.GE.2) THEN 
        CALL DAXPY(L2,-one,WRKESP,1,DA,1)
        DO I = 1, L2
          DA(I) = -DA(I)
        END DO
      ELSE
        CALL DCOPY(L2,WRKESP,1,DA,1)
      END IF

      RETURN
      END
C*MODULE FMOESD  *DECK ESLAGRN
      SUBROUTINE ESLAGRN(l1,l1i,l1j,l1k,nai,naj,nak,NATI,NATJ,NATK,
     *                   orbxch,iodexch,jodexch,kodexch,wrk,wrk1,
     *                  layfrg,idmrec0,ESPA,mapi,mapj,mapk,irec0,iaglob,
     *                   fmode)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      Parameter (half=0.5D+00,zero=0.0D+00)
      logical GOPARR,DSKWRK,MASWRK,doddcor
      logical isini, isinj, isink, orbxch,iodexch,jodexch,kodexch
CZ
      dimension layfrg(*),idmrec0(*)
      DIMENSION ESPA(*),WRK(*),WRK1(*)
      DIMENSION IAGLOB(*)
      DIMENSION mapi(*),mapj(*),mapk(*) 
      DIMENSION FMODE(3,*)

      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      common /fmopnt/ lichfg,lmulfg,lidmrec,lfrgnam,llayfrg,lindat,
     *                lnCBS,lfmozan,lfmoc,lfmomas,lizbas,liaglob,libdgh,
     *                liabdfg,ljabdfg,lnCAO,lidxCAO,liaprjo,ljaprjo,
     *                lCoreAO,lOccCor,lshiftb,liodfmo,lfmoda,lfmodb,
     *                lfmoespa,lfmoespb,llocfmo,lscffrg,lrij,lpopmul,
     *                lpopmat,lialoc,lindbd,liatfrg,lindfrg,lindgfrg,
     *                lnatfrg,lnat0frg,lianfrg,lzanfrg,lcfrg,llibish,
     *                llibnsh,llibng,lindatg,lfmobuf(3),lfmode,lnumfrg,
     *                lloctat,liaoglob,lloadm,lfmoge,ldgrid,liodcfmo,
     *                ljob2grp,lfmopg,lemocdr,luntxyz,luntrot,lstonep,
     *                lmapsu,lfrgmul
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
      PARAMETER (TWO=2.0D+00)
CZ
CZ    THIS CALCULATES LAGRANGIAN-LIKE TERMS (DI+DJ)*VIJ
CZ
      IF (IFMOSTP.NE.4) THEN
        WRITE(IW,*) 'ESLAGRN MUST BE CALLED IN ESD'
        CALL ABRT
      END IF

      ifg = icurfg
      jfg = jcurfg
      kfg = kcurfg
      if (ifg.eq.0.or.jfg.eq.0) then
        write(6,*) 'detect a bug in dimer ESLAGRN'
        call abrt
      end if 
CZ
CZ    L1 is supposed to be for dimers
CZ
      L1A = L1
      doddcor=iand(ixesp,1024).ne.0
CZ
CZ    CALCULATION FOR IFG
CZ
      ilay=icurlay
      klay=min(ilay,layfrg(ifg))
      irec0 = idmrec0(klay)
      if(doddcor) then
        irec0rhf=1
        if(irec0.eq.1) irec0rhf=nfg*2+1
        irec0=irec0rhf
      endif
      idmrec0i=ifg+irec0
CZCZ  call CLOSDA('DELETE')
CZCZ  CALL OPENDA(0)
CZCZ  call makemol(ifg,0,0,klay,0,0,0,0,0,0,0,.true.)
      call makemol(ifg,0,0,klay,0,0,0,0,0,0,0,.false.)
c
      l2i = (l1i*l1i + l1i)/2
      l2j = (l1j*l1j + l1j)/2
      MAXINT = MAX(3*NATI*L2I,3*NATJ*L2J)
      IF (IFMOSTP.EQ.8) THEN
        l2k = (l1k*l1k + l1k)/2
        MAXINT = MAX(MAXINT,3*NATK*L2K)
      END IF
      CALL VALFM(LOADFM)
      LDSINT = LOADFM + 1
      LAST   = LDSINT + MAXINT
      NEED   = LAST - LOADFM - 1
      CALL GETFM(NEED)
CZ 
CZ    HERE THE FIRST ARGUMENT IS DUMMY
CZ
      ncurs  = ncursh
      ncursh = 0
      CALL EXTSDER(WRK,X(LDSINT),L1I,L2I,.TRUE.,.TRUE.,0)
      ncursh = ncurs
      IF(GOPARR) CALL DDI_GSUMF(1501,X(LDSINT),L2I*3*NATI)
c     if (maswrk) write(6,'(8F12.8)') (X(LDSINTI+L2I+I-1),I=1,L2I)

c     call rareads(IDAFMO,x(liodfmo),WRK,L2I,idmrec0i,0)
      call readmond(wrk,orbxch,iodexch,nai,l1i,x(liodfmo),idmrec0i)
c     write(6,*) 'wwwdens',ifg,irec0,WRK(1),wrk(2)

      CALL CPYTSQ(WRK,WRK(L2I+1),L1I,1)
      IC1 = 0
      DO I = 1, NATI
        IATG=IAGLOB(I)
CZ      if (maswrk) write(6,*) 'IATG', IATG
        DO K = 1, 3
          CALL TFTRI(WRK1,X(LDSINT+IC1*L2I),WRK(L2I+1),
     *               WRK1(L2I+1),L1I,L1I,L1I)
          CALL DSCAL(L2I, HALF,WRK1,1)
          IC1 = IC1 + 1

          decon = zero
          loop = 0
          do ii = 1, L1A
            do jj = 1, ii
              temp = zero
              loop = loop + 1
              ii1=mapi(ii)
              ij1=mapi(jj)
              ii2=max(ii1,ij1)
              ij2=min(ii1,ij1)
              isini=ii1.ne.0.and.ij1.ne.0
              if (isini) then
                temp = WRK1((ii2*ii2-ii2)/2+ij2)
              end if

              if (temp.ne.zero) then
                if (ii.eq.jj) then
                  decon = decon +     temp*ESPA(loop)
                else
                  decon = decon + two*temp*ESPA(loop)
                end if
              end if
            end do
          end do

          if (maswrk) FMODE(K,IATG) = FMODE(K,IATG) + decon
        END DO
      END DO

CZ
CZ    CALCULATION FOR JFG
CZ
      klay=min(ilay,layfrg(jfg))
      irec0 = idmrec0(klay)
      if(doddcor) then
        irec0rhf=1
        if(irec0.eq.1) irec0rhf=nfg*2+1
        irec0=irec0rhf
      endif
      idmrec0j=jfg+irec0
CZCZ  call CLOSDA('DELETE')
CZCZ  CALL OPENDA(0)
CZCZ  call makemol(jfg,0,0,klay,0,0,0,0,0,0,0,.true.)
      call makemol(jfg,0,0,klay,0,0,0,0,0,0,0,.false.)

      ncurs  = ncursh
      ncursh = 0
      CALL EXTSDER(WRK,X(LDSINT),L1J,L2J,.TRUE.,.TRUE.,0)
      ncursh = ncurs
      IF(GOPARR) CALL DDI_GSUMF(1501,X(LDSINT),L2J*3*NATJ)
c     if (maswrk) write(6,'(8F12.8)') (X(LDSINT+L2J+I-1),I=1,L2J)

c     call rareads(IDAFMO,x(liodfmo),WRK,L2J,idmrec0j,0)
      call readmond(wrk,orbxch,jodexch,naj,l1j,x(liodfmo),idmrec0j)

      CALL CPYTSQ(WRK,WRK(L2J+1),L1J,1)
      IC1 = 0
      DO J = 1, NATJ
        JATG=IAGLOB(J)
        DO K = 1, 3
          CALL TFTRI(WRK1,X(LDSINT+IC1*L2J),WRK(L2J+1),
     *               WRK1(L2J+1),L1J,L1J,L1J)
          CALL DSCAL(L2J, HALF,WRK1,1)
          IC1 = IC1 + 1

          decon = zero
          loop  = 0
          do ii = 1, L1A
            do jj = 1, ii
              temp = zero
              loop = loop + 1
              ji1=mapj(ii)
              jj1=mapj(jj)
              ji2=max(ji1,jj1)
              jj2=min(ji1,jj1)
              isinj=ji1.ne.0.and.jj1.ne.0
              if (isinj) then
                temp = WRK1((ji2*ji2-ji2)/2+jj2)
              end if

              if (temp.ne.zero) then
                if (ii.eq.jj) then
                  decon = decon +     temp*ESPA(loop)
                else
                  decon = decon + two*temp*ESPA(loop)
                end if
              end if
            end do
          end do

          if (maswrk) FMODE(K,JATG) = FMODE(K,JATG) + decon
        END DO
      END DO

CZ
CZ    CALCULATION FOR KFG
CZ
      IF (IFMOSTP.EQ.8) THEN
        klay=min(ilay,layfrg(kfg))
        irec0 = idmrec0(klay)
        if(doddcor) then
          irec0rhf=1
          if(irec0.eq.1) irec0rhf=nfg*2+1
          irec0=irec0rhf
        endif
        idmrec0k=kfg+irec0
CZCZ    call CLOSDA('DELETE')
CZCZ    CALL OPENDA(0)
CZCZ    call makemol(kfg,0,0,klay,0,0,0,0,0,0,0,.true.)
        call makemol(kfg,0,0,klay,0,0,0,0,0,0,0,.false.)

        ncurs  = ncursh
        ncursh = 0
        CALL EXTSDER(WRK,X(LDSINT),L1K,L2K,.TRUE.,.TRUE.,0)
        ncursh = ncurs
        IF(GOPARR) CALL DDI_GSUMF(1501,X(LDSINT),L2K*3*NATK)
c       if (maswrk) write(6,'(8F12.8)') (X(LDSINT+L2K+I-1),I=1,L2J)

c       call rareads(IDAFMO,x(liodfmo),WRK,L2K,idmrec0k,0)
        call readmond(wrk,orbxch,kodexch,nak,l1k,x(liodfmo),idmrec0k)

        CALL CPYTSQ(WRK,WRK(L2K+1),L1K,1)
        IC1 = 0
        DO LL = 1, NATK
          KATG=IAGLOB(LL)
          DO K = 1, 3
            CALL TFTRI(WRK1,X(LDSINT+IC1*L2K),WRK(L2K+1),
     *                 WRK1(L2K+1),L1K,L1K,L1K)
            CALL DSCAL(L2J, HALF,WRK1,1)
            IC1 = IC1 + 1

            decon = zero
            loop  = 0
            do ii = 1, L1A
              do jj = 1, ii
                temp = zero
                loop = loop + 1
                ki1=mapk(ii)
                kj1=mapk(jj)
                ki2=max(ki1,kj1)
                kj2=min(ki1,kj1)
                isink=ki1.ne.0.and.kj1.ne.0
                if (isink) then
                  temp = WRK1((ki2*ki2-ki2)/2+kj2)
                end if

                if (temp.ne.zero) then
                  if (ii.eq.jj) then
                    decon = decon +     temp*ESPA(loop)
                  else
                    decon = decon + two*temp*ESPA(loop)
                  end if
                end if
              end do
            end do

          if (maswrk) FMODE(K,KATG) = FMODE(K,KATG) + decon
        END DO
      END DO
 
      END IF
      CALL RETFM(NEED)

CZ
CZ    RESTORE THE DIMER CALCULATION
CZ
CZCZ  call CLOSDA('DELETE')
CZCZ  CALL OPENDA(0)
      call makemol(ifg,jfg,kfg,ilay,0,0,0,0,0,0,0,.true.)
      RETURN
      END
