C 14 Nov 05 - DGF - various changes for FMO 2.1 release
C 19 Sep 05 - MWS - kill h or i function runs
C  5 Jul 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  1 Jun 05 - DGF - fixes for the 2nd release
c 15 Mar 05 - DGF - major changes for the second release
C 13 feb 05 - MWS - pad common block NSHEL
C  5 Feb 05 - DGF - FMOESP2,FMO2EI: adaptation to new integral codes
C 10 Nov 04 - MWS - ESD2DER: no need to use GAMGEN anymore
C  7 Sep 04 - MWS - FMO2EI,ESD2DER: tweak rys/rotated axis int. choice
C 19 May 04 - DGF,KK - implement Fragment Molecular Orbital (FMO) method
c
C*MODULE fmoint  *DECK fmosavb
      SUBROUTINE fmosavb(modef,libish,libnsh,libng)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*8 ATOMNM
c     character*4 atm
c     logical MMONLY,QMMM
      PARAMETER (MXSH=5000,MXGTOT=20000,MXATM=2000,MXAO=8192,MAXL=5,
     *           MAXNZ=137,MaxLay=5,
     *           MXSFMO=MaxLay*40,MXGFMO=MaxLay*100,MXAFMO=MaxLay*10)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),ch(mxgtot),ci(mxgtot),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
c     COMMON /TINOPT/ MMONLY, QMMM
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmoshl/ natl,numl,fzan(MXAFMO),fC(MXGFMO,MAXL),fEX(MXGFMO)
     *               ,lshell,LSTART(MXSFMO),LATOM(MXSFMO),LTYPE(MXSFMO),
     *                LNG(MXSFMO),LMIN(MXSFMO),LMAX(MXSFMO),llay(MXAFMO)
      dimension libish(MAXNZ,maxbas,*),libnsh(MAXNZ,maxbas,*),
     *          libng(MAXNZ,maxbas,*)
      equivalence (ATOMNM,FATOMNM)
c
c     save the shell data
c     modef tells where to take the shell data:
c           0: from COMMON /NSHEL / (copy NSHEL to fmoshl)
c       .ne.0: from /fmoshl/ (so that only indexing arrays are filled).
c
      if(modef.eq.0) then
      if(NSHELL.gt.MXSFMO.or.nat.gt.MXAFMO) then
        write(6,*) 'Internal buffer too small',NSHELL,MXSFMO,nat,MXAFMO
        call abrt
      endif
      NCOEF=MXGFMO
cnba+ this is a bug for large basis sets. One has to check if MXGFMO
c     contains all exponents in EX.
c
      CALL DCOPY(NCOEF,CS,1,fC(1,1),1)
      CALL DCOPY(NCOEF,CP,1,fC(1,2),1)
      CALL DCOPY(NCOEF,CD,1,fC(1,3),1)
      CALL DCOPY(NCOEF,CF,1,fC(1,4),1)
      CALL DCOPY(NCOEF,CG,1,fC(1,5),1)
      CALL DCOPY(NCOEF,EX,1,fex,1)
c
      CALL DCOPY(nat,zan,1,fzan,1)
C
      NK=NSHELL
C
      CALL ICOPY(NK,KSTART,1,lSTART,1)
      CALL ICOPY(NK,KATOM ,1,lATOM,1)
      CALL ICOPY(NK,KTYPE ,1,lTYPE,1)
      CALL ICOPY(NK,KNG   ,1,lNG,1)
c     CALL ICOPY(NK,KLOC  ,1,lLOC,1)
      CALL ICOPY(NK,KMIN  ,1,lMIN,1)
      CALL ICOPY(NK,KMAX  ,1,lMAX,1)
c
      natl=nat
      numl=num
      lshell=nshell
c
      do ii=1,nat
        FATOMNM=ANAM(ii)
        iposm=ifndchr(ATOMNM,8,'-')
        iposd=ifndchr(ATOMNM,8,'.')
c       ipos=3
c       if(ATOMNM(2:2).eq.'-') ipos=2
c       atomnm4=ATOMNM(1:ipos-1)
c       iatz=izsymnum(atomnm4)
        ipose=8
        if(iposm.ne.0) then
          read(UNIT=ATOMNM(iposm+1:ipose),FMT='(I3)') ilayer
          if(ilayer.gt.nlayer.or.ilayer.le.0) then
           write(6,*) 'Invalid layer in the basis library',ilayer,nlayer
            call abrt
          endif 
          ipose=iposm-1
        else 
          ilayer=1
        endif 
        if(iposd.ne.0) then
          read(UNIT=ATOMNM(iposd+1:ipose),FMT='(I3)') ibas
          if(ibas.le.0.or.ibas.gt.maxbas) then
            write(6,*) 'Invalid basis number in the basis library',ibas
            call abrt
          endif
        else
          ibas=1
        endif
        llay(ii)=ilayer+ibas*100
c       ilay looks like:
c       201 second basis for layer 1
      enddo
      endif
c
c     find the maximim number of AOs on a single atom
c     fill indexing arrays, while assuming that each atom's basis set is
c     stored continuously.
c
      call viclr(libish,1,MAXNZ*maxbas*nlayer)
      call viclr(libnsh,1,MAXNZ*maxbas*nlayer)
      call viclr(libng,1,MAXNZ*maxbas*nlayer)
      maxcbs=0
      jat=0
      jbas=0
      jlay=0
      nao=0
c     mxangm=0
      do ii=1,lshell
        iat=latom(ii)
        ilay=mod(llay(iat),100)
        ibas=llay(iat)/100
        iz=int(fzan(iat)+0.1D+00)
        if(iat.ne.jat.or.jlay.ne.ilay.or.jbas.ne.ibas) then
           maxcbs=max(maxcbs,nao)
           nao=0
           jat=iat
           jbas=ibas
           jlay=ilay
           if(libish(iz,ibas,ilay).ne.0) then
             write(iw,9000) iz,ibas,ilay 
             call abrt
           endif
           libish(iz,ibas,ilay)=ii 
        endif
c       mxangm=max(mxangm,KTYPE(ii)-1)
        nao=nao+lmax(ii)-lmin(ii)+1
        libnsh(iz,ibas,ilay)=libnsh(iz,ibas,ilay)+1
        libng(iz,ibas,ilay)=libng(iz,ibas,ilay)+lng(ii)
      enddo
      maxcbs=max(maxcbs,nao)
      return
 9000 format(1x,'Redundant basis sets for iat=',I2,', ibas=',I2,
     *          ', ilay=',I2)
      end
c
C*MODULE fmoint  *DECK fmoesp
      SUBROUTINE fmoesp(l1,l2,h,layfrg,scffrg,idmrec0)
      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,
     *           one=1.0D+00,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
      dimension h(*),layfrg(*),scffrg(*),idmrec0(*),karten(0:maxl-1)
      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)
c     COMMON /MP2PAR/ OSPT,TOL,METHOD,NWDMP2,MPPROP,
c    *                NACORE,NBCORE,NOA,NOB,NORB,NBF,NOMIT
      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)
c     COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      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,convfg,cnvdmp,coroff,IXESP,mxitfg,nguess,
     *                NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c
      common /fmopnt/ lichfg,lmulfg,lidmrec,lfrgnam,llayfrg,lindat,
     *                lnCBS,lfmozan,lfmoc,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(2),lfmode,lnumfrg,lloctat,
     *                liaoglob,lloadm,lfmoge,ldgrid,liodcfmo
      common /fmorun/ espscf,IDAFMO,icurlay,icurfg,jcurfg,kcurfg,ncursh,
     *                ngau,icurpop,ifmostp,moncor,nat1e,needr,modrst,
     *                norbproj
      data karten/1,4,6,10,15/ 
      data rhf/8HRHF     /,RMC/8HMCSCF   /,
     *     dbgfmo/8HDBGFMO  /,dbgme/8HFMOESP  /,debug/8HDEBUG   /
c
c     this subroutine 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
      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
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)
      esdder=ifmostp.eq.6.and.resdim.ne.0.and.iand(ixesp,32).eq.0.and.
     *       nder.gt.0
      jfg=jcurfg
      lfg=kcurfg
      if(ifmostp.eq.6) then
        kfg0=jfg
        kfg1=kfg0
        jfg=0
        iz=2
        kfact=nat1e
      else
        kfg0=1
        kfg1=nfg
        iz=1
        kfact=1
      endif 
      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
      if(iz.eq.2) einuc=ENUC(nat,Zan,C)
c
      call vclr(x(lfmoespa),1,l2)
cb    call vclr(x(lfmoespb),1,l2)
c
      if(iz.eq.1) CALL DAread(IDAF,IODA,H,L2,11,0)
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
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
          if(iz.eq.1.and.rk.gt.resppci.and.resppci.ne.zero) goto 100
          espap=rk.gt.respapi.and.respapi.ne.zero
        else
          espap=.false.
        endif
c
        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 ESP 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).
        if(iz.eq.2) then
c         set nat1e to jfg atoms
          nat1e=nat-nat0 
c         basis set must be reset to ifg monomer to do 1e integrals. 
c         save dimer info
          nshs=nshell
          nums=num
          nats=nat
          num=num0
          nshell=ncursh
          nat=nat0
c         icurfg=0 prevents unwanted recursive calls to FMOESP from ONEEI.
c         icurfg is stored in ncursh, that is unused in oneei.
          ncurs=ncursh
          ncursh=icurfg 
          icurfg=0
c
          call oneei
c
          ncursh=ncurs
          nat1e=natfmo
          num=nums
          nshell=nshs
          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
          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         contributions. One could equally write ifg.lt.jfg.
c
        if(ifg.gt.jfg) epot=ENUCR-einuc-ENUC(nat-nat0,Zan(n0j),C(1,n0j))
          if(esdder.and.kfact.eq.0) then
c           Save 1e potential 
            call dcopy(l2,x(lfmoespb),1,x(lwrkesp),1)
          endif
          if((rk.gt.resppci.and.resppci.ne.zero).or.kfact.eq.0) then
            goto 100
          endif
        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
            call rareads(IDAFMO,x(liodfmo),x(lfmoda),l2k+l3k,idmrec0k,0)
          else
            call rareads(IDAFMO,x(liodfmo),x(lfmoda+l2k),l3k,idmrec0k,0)
            call DMTX2(x(lfmoda),x(lfmoda+l2k),nak,l1k,l1k) 
          endif
        else
          call rareads(IDAFMO,x(liodfmo),x(lfmoda),l2k,idmrec0k,0)
        endif
c       goto 100
c       write(6,*) 'original Density',kfg,l1k,nak
c       call prtri(X(lfmoda),l1k)
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
          ldena=lfmoda 
          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
            CALL SHLDEN(scftyp1,X(lwrkden),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
        if(esdder) call vclr(x(lwrkesp2),1,l2k)
c       call timit(1)
c       do iohno=1,100
        CALL fmo2ei(SCHWRZ,NINT,NSCHWZ,NSCHWZB,NSCHWNZB,L1,L2,X(LXINTS),
     *              NSH2,X(LGHOND),MAXG,IA,ldena,ldenp,X(lfmoespa),
     *              x(lwrkden),x(lwrkesp2),X(LDSH),X(LDSHb),x(liaglob),
     *              x(lindat),x(lindatg),bssedim,respapi,resppci,espap,
     *              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)
        CALL RETFM(NEED)
c
        if(esppar) goparr=savgop
        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)
          CALL daxpy(L2k,one,X(lwrkesp2),1,X(lwrkesp),1)
        endif
  100 continue
      IF(esppar.and.goparr.and.nxt) CALL DDI_DLBRESET
C
C   --- OFF DIAGONAL ELEMENTS ARE DOUBLE THE CORRECT VALUE ---
C
      espscf0=espscf 
      if(ifmostp.eq.6) espscf0=kfact*half
      if(kfact.ne.0) then
        CALL DSCAL(L2,half*espscf0,X(lfmoespa),1)
        II=lfmoespa-1
        DO 220 I=1,L1
          II = II+I
          X(II) = X(II) + X(II) 
  220   CONTINUE
c
c       Sum over nodes the ESP matrix.
c
        if(goparr) call ddi_gsumf(2418,x(lfmoespa),l2)
c       call ddi_gsumf(2419,x(lfmoespb),l2)
      endif
c
c     if(ifmostp.eq.4) some=.true.
      if(some) then
        write(iw,*) 'Final 1e ESP is'
        call prtril(x(lfmoespb),num0)
        write(iw,*) 'Final 2e ESP is'
        call prtril(x(lfmoespa),num0)
      endif
c     Add 1e ESP stored in lfmoespb.
      if(ifmostp.eq.6) espscf0=one
      call daxpy(l2,espscf0,x(lfmoespb),1,x(lfmoespa),1)
      if(some) then
        write(iw,*) 'Final 1e+2e ESP is'
        call prtril(x(lfmoespa),num0)
      endif
      if(esdder.and.kfact.ne.0) then
        call ESDnrder(nat0)
c       write(6,*) 'density I',l1
c       call prtril(x(lwrkden),l1)
c       write(6,*) 'density J',l1k
c       call prtril(x(lfmoda),l1k)
        CALL esd2der(x(lwrkden),x(lfmoda),l1)
c       add the contributions (NR+2e) to the gradient
c       Not only masters accumulate the gradient here: all nodes have
c       contributions (no global sum in the above two calls).
        call fmodeg(1,x(lfmode+natfmo*3),x(liaglob))
      endif
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
c
c       finally, add ESP to the one-electron Hamiltonian
c
        if(iz.eq.1) then
          call daxpy(l2,one,x(lfmoespa),1,h,1)
          CALL DAwrit(IDAF,IODA,h,L2,11,0)
        endif
      endif
      nprfmo=nprsav
      dirscf=dirsav
c     write(6,*) 'Exit of FMOESP',nqmt
      return
 9000 format(/1x,'ESP lfg=',I5,': NZ',I12,' skipped',I10,' blocks.')
 9005 format(/1x,'ESP lfg=',I5,': NZ',I12,' skipped',I10,' blocks',I8,
     *           ' superblocks.')
      end
C*MODULE fmoint  *DECK fmoesp2
      SUBROUTINE fmoesp2(IA,DA,FA,DB,FB,GHONDO,l1,NINT,iloop,espap,
     *                   esdder)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL IANDJ,KANDL,SAME,OUT,espap,esdder,espap3,espap4
C
      DIMENSION IA(*),DA(*),FA(*),DB(*),FB(*),GHONDO(*),inds(4)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, four=4.0D+00)
C
      COMMON /ERIOUT/ ISH,JSH,KSH,LSH,LSTRI,LSTRJ,LSTRK,LSTRL
      COMMON /MISC  / IANDJ,KANDL,SAME
      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 /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJX,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
      common /fmorun/ espscf,IDAFMO,icurlay,icurfg,jcurfg,kcurfg,ncursh,
     *                ngau,icurpop,ifmostp,moncor,nat1e,needr,modrst,
     *                norbproj
C
      DATA HALF /0.5D+00/
c     DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF     /
c     DATA RMC/8HMCSCF   /,GVB/8HGVB     /
C
c     This subroutine computes 2e contributions to the ESP on a given 
c     monomer(dimer) by an external monomer. Density matrices are from
c     that external monomer.
c     In (II JJ | KK LL) II and JJ are from the external monomer;
c                        KK and LL are from the main monomer(dimer). 
c     any of II,JJ are NEVER equal to any of KK,LL.
c     There is much less permutation symmetry compared to the usual case 
c     when all 4 indices are from the same entity. 
c     Usually only Coulomb contribution is added (iloop.eq.0). 
c     Exchange contribution is added with iloop.eq.1.
c     cloned from DIRFCK.
c     
c     parstat: NodeNone 
C
C     NOTE THAT OFF-DIAGONAL ELEMENTS WILL NEED TO BE HALVED LATER.
C
c     UROHF = SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
c     ALFSCF= SCFTYP.EQ.RHF  .or.  scftyp.eq.rmc
      CUTINT = CUTOFF
C
c     IF(SCFTYP.EQ.GVB) THEN
c        NSHL = NHAM
c        IF(NCO.GT.0) NSHL=NSHL-1
c     END IF
C
      SAME = ISH .EQ. KSH .AND. JSH .EQ. LSH
      IANDJ = ISH .EQ. JSH
      KANDL = KSH .EQ. LSH
      MINI = KMIN(ISH)
      MINJ = KMIN(JSH)
      MINK = KMIN(KSH)
      MINL = KMIN(LSH)
      MAXI = KMAX(ISH)
      MAXJ = KMAX(JSH)
      MAXK = KMAX(KSH)
      MAXL = KMAX(LSH)
      LOCI = KLOC(ISH)-MINI
      LOCJ = KLOC(JSH)-MINJ
      LOCK = KLOC(KSH)-MINK
      LOCL = KLOC(LSH)-MINL
c
c     Find which 2 out of 4 indices come from the external monomer. These would
c     normally be ISH and JSH (see fmo2ei), but some integral drivers 
c     (Pople and Eric) like to reorder integrals according to the L value.
c     Such reordering is irrelevant for ordinary integrals in RHF but has to be
c     untangled here.
c
        indi=1
        indj=2
        indk=3
        indl=4
      if(ish.gt.ncursh.and.jsh.gt.ncursh) then
c       set above 
      else if(ish.gt.ncursh.and.ksh.gt.ncursh) then
        indi=1
        indj=3
        indk=2
        indl=4
      else if(ish.gt.ncursh.and.lsh.gt.ncursh) then
        indi=1
        indj=4
        indk=2
        indl=3
      else if(jsh.gt.ncursh.and.ksh.gt.ncursh) then
        indi=2
        indj=3
        indk=1
        indl=4
      else if(jsh.gt.ncursh.and.lsh.gt.ncursh) then
        indi=2
        indj=4
        indk=1
        indl=3
      else if(ksh.gt.ncursh.and.lsh.gt.ncursh) then
        indi=3
        indj=4
        indk=1
        indl=2
      else
c       We do not really check if more than 2 are .gt.ncursh but that is
c       not needed as long as fmo2ei is sane. 
        call abrt
      endif
      espap3=espap.and.indj.eq.3
      espap4=espap.and.indj.eq.4
c
c     espap is a headache fountain.
c     For espap we want to do diagonal terms in the second external monomer
c     shell. Such shell is stored in indj 
c     (that is, indj=2 means JSH, indj=3 is KSH and indj=4 is LSH).
C
      IJN = 0
      jmin=minj
      JMAX = MAXJ
      DO 360 I = MINI,MAXI
         I_INDEX = (I-MINI)*LSTRI + 1
         IF (IANDJ) JMAX = I
         inds(1) = LOCI+I
         if(espap.and.indj.eq.2) then
c          ESPAP does only diagonal terms for the 2nd index (1st=2nd).
           jmin=jmax
           IJN=IJN+jmin-minj 
         endif
         indv1=i
         DO 340 J = jmin,JMAX
            IJN = IJN+1
            IJ_INDEX = (J-MINJ)*LSTRJ + I_INDEX
c           fix index offset to point to the proper elements
            inds(2) = LOCJ+J
            LMAX = MAXL
            KLN = 0
            if(indi.eq.2) indv1=j
            DO 320 K =  MINk,maxk
               IJK_INDEX = (K-MINK)*LSTRK + IJ_INDEX
               IF (KANDL) LMAX = K
               inds(3) = LOCK+K
               if(indi.eq.3) indv1=k
               DO 300 L = minl,LMAX
                  KLN = KLN+1
c                 skip off-diagonal terms for espap.
c                 Clumsy but perhaps fairly clear.
c                 indv1 points to the first index (i,j or k) with which 
c                 the second (k or l) has to be equal. 
c                 The whole fuss is raised to cover the knavish SAME. 
                  if(espap3.and.k.ne.indv1 .or. espap4.and.l.ne.indv1)
     *              goto 300
                  if(SAME .AND. KLN.GT.IJN) GO TO 340
                  IJKL_INDEX = (L-MINL)*LSTRL + IJK_INDEX
C
                  VAL = GHONDO( IJKL_INDEX )
                  IF(ABS(VAL).LT.CUTINT) GO TO 300
                  NINT = NINT + 1
C
                  inds(4) = LOCL+L
                  if(iloop.eq.0) then
                     ii = inds(indi)
                     jj = inds(indj)
                     kk = inds(indk)
                     ll = inds(indl)
c                    ii,jj,kk,ll are now in the semicanonic order. They only
c                    need to be reordered triangularly (I>J).
c                    ii,jj are external, kk,ll are n-mer.
                     if(ii.lt.jj) then
                       n=ii
                       ii=jj
                       jj=n
                     endif
                     if(kk.lt.ll) then
                       n=kk
                       kk=ll
                       ll=n
                     endif
                     if(espap) then
                       NIJ=II-l1
                     else
                       NIJ = IA(II-l1)+JJ-l1
                     endif
                     NKL = IA(KK)+LL
                     IF(II.EQ.JJ) VAL = VAL*HALF
                     IF(KK.EQ.LL) VAL = VAL*HALF
                     IF(II.EQ.KK.AND.JJ.EQ.LL) VAL = VAL*HALF
                     VAL4 = VAL*four
                  else
c
c                    to the best of the present knowledge the code below is
c                    correct (even for Eric). However, if ixesp=1 is to be
c                    used for something real, perhaps detailed comparison
c                    (i.e. RHF exchange vs exchange here) should be conducted.
c
                     ii = inds(indi)
                     jj = inds(indk)
                     kk = inds(indj)
                     ll = inds(indl)
c                    ii,jj,kk,ll are now in the semicanonic order. They only 
c                    need to be reordered triangularly (I>J).
c                    ii,kk are external, jj,ll are n-mer.
                     if(ii.lt.kk) then
                       n=ii
                       ii=kk
                       kk=n
                     endif
                     if(jj.lt.ll) then
                       n=jj
                       jj=ll
                       ll=n
                     endif
                     if(espap) then
                       NIJ=II-l1
                     else
                       NIJ = IA(II-l1)+kk-l1
                     endif
                     NKL = IA(jj)+LL
                     IF(II.EQ.JJ) VAL = VAL*HALF
                     IF(KK.EQ.LL) VAL = VAL*HALF
                     IF(II.EQ.KK  .AND.  JJ.EQ.LL) VAL = VAL*HALF
                     VAL4 = -VAL
c                    Exchange is -1/2 compared to Coulomb.
                  endif
C
C      WE NOW HAVE EVERYTHING READY, PUT INTEGRALS INTO FOCK MATRIX
C
c                 external monomer had shell indices in the lower part
c                 of the basis set storage, density matrix, however, is
c                 stored as if the monomer were the only thing on Earth.
C
                  FA(NKL) = FA(NKL)+VAL4*DA(NIJ)
                  if(esdder) then  
                     FB(NIJ) = FB(NIJ)+VAL4*DB(NKL)
                  endif 
  300          CONTINUE
  320       CONTINUE
  340    CONTINUE
  360 CONTINUE
      RETURN
      END
C*MODULE fmoint  *DECK fmo2ei
      SUBROUTINE fmo2ei(SCHWRZ,NINT,NSCHWZ,NSCHWZB,NSCHWNZB,L1,L2,XINTS,
     *                  NSH2,GHONDO,MAXG,IA,lDA,ldp,FA,DB,FB,DSH,DSHb,
     *                  iaglob,indat,indatg,bssedim,respapi,resppci,
     *                  espap,esdder,ifg,jfg,kfg,lfg)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL OUT,SCHWRZ,GOPARR,DSKWRK,MASWRK,NXT,bimer(3),bssedim,
     *        espap,espapij,smartr(2),esdder
C
      DIMENSION XINTS(NSH2),GHONDO(MAXG),IA(L1),FA(L2),DB(*),FB(*),
     *          DSH(NSH2),DSHb(NSH2),iaglob(*),indat(*),indatg(natfmo,*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
c     COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
c    *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),ch(mxgtot),ci(mxgtot),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
c     COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
c     COMMON /PKFIL / PK,PANDK,BLOCK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      common /shlexc/ norgsh(3),norgsp(3),iexch,nangm,ngth(4)
      COMMON /SHLG70/ ISH,JSH,KSH,LSH,IJKLXX(4)
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
      COMMON /FMCOM / X(1)
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,convfg,cnvdmp,coroff,IXESP,mxitfg,nguess,
     *                NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c
      common /fmorun/ espscf,IDAFMO,icurlay,icurfg,jcurfg,kcurfg,ncursh,
     *                ngau,icurpop,ifmostp,moncor,nat1e,needr,modrst,
     *                norbproj
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
c     integer ida(2) 
C
      PARAMETER (ZERO=0.0D+00,FOUR=4.0D+00)
C
      data bimer/3*.false./
C
C     Compute 2e contribution to ESP for FMO. This subroutine contains the
c     outer 4 shell loops and it is cloned from TWOEI. Additional comments are
c     available there and most of them apply.  Much less symmetry 
c     (no 1-2 particle permutation symmetry; enforced C1 symmetry)
c     promted not calling TWOEI but having an independent subroutine. 
c     For smart distances one- and two-electron terms must be used in the
c     same fashion, based on these smart distances!
c     PK option is NOT supported!
c     parstat: GroupNone/GroupFull 
C
      TIM = ZERO
      CALL TSECND(TIM)
C
C     ----- INITIALIZATION FOR PARALLEL WORK -----
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT = -1
      MINE = -1
c     CMBDIR= DIRSCF .OR. DIRNLO .OR. DIRTRF
C
      ngth(4) = 1 
      ngth(3) = ngth(4) * NANGM
      ngth(2) = ngth(3) * NANGM
      ngth(1) = ngth(2) * NANGM
      do i=1,3
         norgsh(i) = 0
         norgsp(i) = 0
      enddo
C
      NINT  = 0
      NSCHWZ= 0
      NSCHWZB= 0
      NSCHWNZB= 0
      DENMAX = ZERO
      smartr(1)=iand(ixesp,8).ne.0.and.jfg.ne.0
      smartr(2)=iand(ixesp,16).ne.0.and.jfg.ne.0
c     turn off smart distances for connected n-mers.
      if(smartr(1).and.nbdfg.ne.0) then
        bimer(1)=fmodist(ifg,0,0,jfg).eq.0
        if(kfg.eq.0) then
          if(bimer(1)) smartr(1)=.false.
        else
          bimer(2)=fmodist(ifg,0,0,kfg).eq.0
          bimer(3)=fmodist(jfg,0,0,kfg).eq.0
        if(bimer(1).and.(bimer(2).or.bimer(3)).or.bimer(2).and.bimer(3))
     *    smartr(1)=.false.
        endif
      endif
      if(schwrz.and..not.esdder) xintmax=xints(idamax(nsh2,xints,1))
c
c     The problem with skipping superblocks of shells and esdder is
c     not fundamental but rather foolish ignorance of trying to do
c     both contractions over 1st and 2nd electrons in one 4-loop,
c     rather than contracting them separately (in which case one could
c     skip superblocks for esdder). Since normally esdder is true
c     (for runs with gradient, that is), this addition of skipping superblocks 
c     is fairly useless, except mostly 3-body runs for which esdder is always 
c     false (with the present code). 
c
c     for FMO ESP runs the indices are divided as follows:
c     particle one (II,JJ): interacting monomer K   (ncursh+1,...,nshell) 
c     particle two (KK,LL): original monomer(dimer) (1,...,ncursh)
c     ncursh gives the border line between the two in the stored basis.
c     K is looped over somewhere above.
c
C     ----- I SHELL -----
C
      jj0=ncursh+1
      DO 920 II = ncursh+1,NSHELL
C
C       ----- CHECK CPU TIME -----
C
        CALL TSECND(TIM)
        IF(TIM.GE.TIMLIM) THEN
          IF(MASWRK) WRITE(IW,9030)
          CALL ABRT
        END IF
C
C       ----- PRINT INTERMEDIATE RESTART DATA -----
C
c       IF(NPRINT.NE.-5  .AND.  .NOT.CMBDIR .AND. MASWRK)
c    *     WRITE(IW,9010) II,JST,KST,LST,NREC,ICOUNT
C
C       ----- J SHELL -----
C
        if(espap) jj0=II 
        DO 900 JJ = jj0,II
C
          IF(SCHWRZ) THEN
            IJIJ = (II*II-II)/2 + JJ
c           DSH is written for the external monomer with no offset
c           espap and not espapij should be used because DSH was written
c           with espap ere shells came to be.
            if(espap) then
              dmaxij=FOUR*dsh(II-ncursh)
            else
              dmaxij=FOUR*dsh(IA(II-ncursh)+JJ-ncursh)
            endif
            if(schwrz.and..not.esdder) then
c             write(6,*) ii,jj,'www2ei',XINTS(IJIJ),xintmax,dmaxij
              if(XINTS(IJIJ)*xintmax*dmaxij.LT.CUTOFF) then
                NSCHWZB=NSCHWZB+1
                goto 900
              endif
            endif
          END IF
c         for smart distances Schwarz tests are based upon full density,
c         not populations if espap is false, even though it is populations that
c         are used below (espapij true). The reason is DSH that has been 
c         constructed before it was known if one uses populations or density.
c         There may be some room for improvement here.
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 900
          END IF
          NSCHWNZB=NSCHWNZB+1
C
C         ----- K SHELL -----
C
          DO 880 KK = 1,ncursh
c
c           write(6,*) 'K checking shel',kk,KATOM(kk),indat(iaglob(KATOM(kk)))
c           if(bssedim.and.indat(iaglob(KATOM(kk))).eq.jcurfg) goto 880
            kkat=iaglob(KATOM(kk))
            kkfg=indat(kkat)
            if(bssedim.and.kkfg.ne.icurfg) goto 880
c           write(6,*) '  K allowed'
C
C           ----- L SHELL ----
C
            DO 860 LL = 1,KK
c
c             write(6,*) 'L checking shel',ll,KATOM(ll),indat(iaglob(KATOM(ll)))
c             if(bssedim.and.indat(iaglob(KATOM(ll))).eq.jcurfg) goto 860
              llat=iaglob(KATOM(ll))
              llfg=indat(llat)
              if(bssedim.and.llfg.ne.icurfg) goto 860
                  if(smartr(1).or.smartr(2)) then
              if(smartr(1)) then
                rk=fmossdist(kkfg,llfg,indatg(kkat,1),indatg(llat,1),
     *                   ifg,jfg,kfg,lfg,bimer)
              else
                rk=fmosdist(kkfg,llfg,indatg(kkat,1),indatg(llat,1),
     *                  ifg,jfg,kfg,lfg)
              endif
              if(rk.gt.resppci.and.resppci.ne.zero) goto 860
c             skip fragments treated by point charges
              espapij=rk.gt.respapi.and.respapi.ne.zero
              if(espap.and..not.espapij) call abrt 
c             double check: espapij is based on larger or equal distance that
c             espap and thus espapij is more often true than espap
c             (approximations are more efficient with smart distances).
              if(espapij.and.ii.ne.jj) goto 860 
c               skip off-diagonal for fragments treated by atomic populations 
c               it is nay good that ii,jj are outer loops, not kk,ll.
              else
                espapij=espap
              endif
              ldd=lda
              if(espapij) ldd=ldp
c             write(6,*) '  L allowed'
C  
C             ----- GO PARALLEL! -----
C  
              IF ((.NOT.NXT) .AND. GOPARR) THEN
               IPCOUNT = IPCOUNT + 1
               IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 860
              END IF
C  
C             ----- (II,JJ//KK,LL) -----
C  
              IEXCH = 1
              ISH = II
              JSH = JJ
              KSH = KK
              LSH = LL
              QQ4 = 1
C  
C             ----- COMPUTE TWO-ELECTRON INTEGRALS ----
C  
C             APPLY THE SCHWARZ INEQUALITY, WHICH IS
C             (II,JJ//KK,LL) .LE.  SQRT( (II,JJ//II,JJ)*(KK,LL//KK,LL) )
C             SEE, FOR EXAMPLE, J.L.WHITTEN, J.CHEM.PHYS. 58,4496-4501(1973)
C      
              IF(SCHWRZ) THEN
                KLKL = (KSH*KSH-KSH)/2 + LSH
                if(esdder) then
                  denmax=max(dmaxij,FOUR*dshb(IA(KSH)+LSH))
                else
                  denmax=dmaxij
                endif
                TEST = XINTS(IJIJ)*XINTS(KLKL)*DENMAX
                if(TEST.LT.CUTOFF) then
                  NSCHWZ = NSCHWZ + 1
                  GOTO 840
                endif
              END IF
c  
              nloop=0
              if(iand(IXESP,1).ne.0) nloop=2
              do 830 iloop=0,nloop
                if(iloop.eq.1) then
c                 compute exchange terms
                  iexch=2
                  ISH = II
                  JSH = kk
                  KSH = jj
                  LSH = LL
                endif
                if(iloop.eq.2) then
                  if(ii.eq.jj.or.ll.eq.kk) goto 830
c                 compute exchange terms
                  iexch=3
                  ISH = II
                  JSH = ll
                  KSH = jj
                  LSH = kk
                endif
c   
                call shellquart(ish,jsh,ksh,lsh,ghondo)
C   
                call fmoesp2(IA,x(lDd),FA,DB,FB,GHONDO,l1,NINT,iloop,
     *                       espapij,esdder)
  830         continue 
C
  840         CONTINUE
  860       CONTINUE
  880     CONTINUE
  900   CONTINUE
  920 CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
      IF(NXT  .AND.  GOPARR) CALL DDI_DLBRESET
      IF(SCHWRZ.and.GOPARR) THEN
         idum=0
         call ddi_nsumi(1055,NSCHWZ,NINT,idum,idum,2)
c        IF(NPRINT.NE.-5 .AND. MASWRK) WRITE(IW,9020) NSCHWZ
      END IF
      RETURN
C
 9030 FORMAT(//1X,'*** THIS JOB HAS EXHAUSTED ITS CPU TIME ***'/
     *         1X,'     (WHILE COMPUTING 2E- INTEGRALS)'///)
      END
C*MODULE fmoint  *DECK fmohop
      subroutine fmohop(l1,l2,h,s,ss,q,dd,scr,wrk1,wrk2,iabdfg,jabdfg,
     *                  idxCAO,iaglob,nCBS,nCAO,iaprjo,japrjo,shiftb,
     *                  CoreAO,fmoc,rotlcao,locfmo,NSHELL,KATOM,KTYPE,
     *                  KLOC,kmin,saveh)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      integer rightend
      logical saveh,locsav,some,GOPARR,DSKWRK,MASWRK,nxt,parrgo,modQAbas
     *       ,modQMbas
      PARAMETER (MXSH=5000,MXATM=2000,one=1.0D+00)
      dimension h(*),s(*),ss(l1,l1),q(l1,l1),dd(*),scr(l1,8),wrk1(*),
     *          wrk2(*),iabdfg(*),jabdfg(*),idxCAO(MaxBnd,*),iaglob(*),
     *          nCBS(*),nCAO(*),iaprjo(MaxCAO,*),japrjo(MaxCAO,*),
     *          CoreAO(MaxCBS,MaxCAO,*),shiftb(MaxCAO,*),fmoc(3,*),
     *          rotlcao(MaxCBS,*),KATOM(*),KTYPE(*),KLOC(*),kmin(*),
     *          locfmo(2,2,*)
      dimension zaxis(3),bond(3)
      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 /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TT(432),INVT(48),NT
      COMMON /SYMREP/ IRPNAM(14),IPA(14),LAMBDA(14),LAMBD0(14),
     *                IADDR1(14),IADDR2(14),IADDR3(14)
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,nsafmo
      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,convfg,cnvdmp,coroff,IXESP,mxitfg,nguess,
     *                NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c
      common /fmorun/ espscf,IDAFMO,icurlay,icurfg,jcurfg,kcurfg,ncursh,
     *                ngau,icurpop,ifmostp,moncor,nat1e,needr,modrst,
     *                norbproj
      data dbgfmo/8HDBGFMO  /,dbgme/8HFMOHOP  /,debug/8HDEBUG   /
      data zaxis/0,0,1/
c
c     Compute hybrid orbital (HO) projector terms that assign a part
c     of usually 5 HO orbitals (1s core and sp3) to a given fragment.
c     iaotyp (same as KTYP in NSHEL) is obsolete now?
c     parstat: GroupFull (to be improved?) 
c
      nsafmo=nsalc
c     nsafmo will store original NSALC before eliminating unwanted
c     AOs for FMO.
      norbproj=0
      if(nbdfg.eq.0) return
      some=(exetyp.eq.debug.or.exetyp.eq.dbgfmo.or.exetyp.eq.dbgme).and.
     *     maswrk
c
      ilay=icurlay
      ihnew=0
      modQAbas=iand(ixesp,128).ne.0.and.saveh
      modQMbas=iand(ixesp,256).ne.0.and.saveh
      parrgo=goparr.and.iand(modpar,16).eq.0.and..not.modQAbas
      locsav=(ifmostp.eq.1.and.irststp.ge.2)
c     For modQbas the metric (SALC etc) is not modified for Huckel
c     (Otherwise some mismatch occurs in COPROJ) and
c     no need to do anything for ESD energy+gradient (ifmostp equal to 6,7).
      if((modQAbas.or.modQMbas).and.
     *   (ifmostp.eq.1.or.ifmostp.eq.6.or.ifmostp.eq.7)) locsav=.true.
c     locsav defines a bizarre run when only locfmo array is filled
c     (as used by some restart jobs).
c     check that no offending code (touching dummies) is put where locsav runs.
      l3=l1*l1
      if(.not.locsav) then
        if(modQAbas) then
c         ss now stores Z, which is 6d -> 5d SALC matrix.
          call daread(IDAF,IODA,ss,L3,44,0)
c         wrk2 stores linear (l1) index for q
          call indsalc(ss,wrk2,nsalc,l1)
        else
          if(saveh) then
            CALL DAread(IDAF,IODA,S,L2,12,0)
            call CPYTSQ(s,ss,l1,1)
            if(modQMbas) then
c             call vclr(h,1,l2)
              write(6,*) 'Orig Hams'
              call prtril(h,l1)
            else
              if(parrgo) then
                CALL DAread(IDAF,IODA,H,L2,11,0)
c               to avoid double counting in parallel have all slaves clear H 
                if(.not.maswrk) call vclr(h,1,l2)
c               ironically DAREAD broadcasts unneeded h to slaves. 
c               Should we fake sequential run by fiddling with GOPARR?
              endif
            endif
          else
c           call RUNITV(l1,l1,ss)
            call vclr(ss,1,l1*l1)
          endif
        endif
      endif
c
C     ----- INITIALIZATION FOR PARALLEL WORK -----
C
      NXT = IBTYP.EQ.1
      NEXT  = -1
      kount = -1
      do 300 ibdg=1,nbdfg
        if(parrgo) then
          kount=kount+1
          IF(NXT) THEN
            IF(kount.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
            if(NEXT.ne.kount) goto 300
          else
            if(MOD(kount,NPROC).NE.me) goto 300
          endif
        endif
c       atoms between which the bond is cut. 
c       ibda=iabdfg(ibdg) 
c       jbda=jabdfg(ibdg) 
c       ibdabs=abs(ibda) 
c       jbdabs=abs(jbda) 
        ierr=0
        if(iabdfg(ibdg).lt.0) then
          leftend=-iabdfg(ibdg)
          rightend=jabdfg(ibdg)
          if(rightend.lt.0) ierr=1 
        else if(jabdfg(ibdg).lt.0) then
          leftend=-jabdfg(ibdg)
          rightend=iabdfg(ibdg)
        else
          ierr=1
        endif
        if(ierr.ne.0) then
          write(iw,*) 'Confusion in FMOHOP:',iabdfg(ibdg),jabdfg(ibdg)
          call abrt
        endif
        ibdtyp=idxcao(ibdg,ilay)
c       find the negative side; that is where the basis set for the overlaps
c       (positive one has a ghost atom added at the negative side so we 
c       always want the negative side).
c       is located.
        ial0=0
        ial=0
        iar=0
c       do not accept bonds that are fully inside a dimer (in which case
c       atomic charges are not modified), that is, only work with split atoms
c       having Z-1 and 1 charges.
        do iat=1,nat 
          izat=int(zan(iat)+1.0D-02)
          if(iaglob(iat).eq.leftend) ial0=iat
          if(izat.ne.ian(iat)) then
            if(iaglob(iat).eq.leftend) ial=iat
            if(iaglob(iat).eq.rightend.and.izat.ne.1) iar=iat
c           izat.ne.1 guards against false propagation when a ghost atom
c           attracts a second broken bond it is involved in.
          endif
        enddo
c       enforce precedence of the left end if both are there and the left end
c       is Z-1. This is neccessary for complicated cases when an atom is
c       involved into two bonds with different ends, e.g.
c       -1 2
c       -2 3
        jat=0
        if(ial.ne.0) then
          jat=ial
        else
          jat=iar
        endif 
c       the left end should be a ghost atom if the right end is in 
c       otherwise the whole bond is inside and it does contribute
c       write(6,*) 'wwwhuhu',jat,iar,ial0,zan(ial0),iat 
        if(iar.ne.0.and.ial0.ne.0.and.int(zan(ial0)+1.0D-02).ne.1) jat=0
        if(jat.ne.0) then
          iside=0
          if(int(zan(jat)+1.0D-02).eq.1) iside=1
          if(some.and..not.locsav) write(6,*) 'Found bond',ibdg,jat
c
c         now find the location where the basis set for the projection 
c         orbitals starts in the overlap matrix
          jj=0
          do ii=1,nshell
            iat=kATOM(ii)
            if(iat.eq.jat) then
              jj=ii
              goto 200
            endif 
          enddo 
          if(maswrk) write(iw,*) 'Bond atom not found',jat,ibdg
          call abrt
  200     continue
          iloc=kloc(jj)
c         save iloc of monomers during monomer initial guess 
c         for dimer density construction
c         The structure of locfmo:
c         locfmo(1,1,ibdg) contains the left side fragment for ibdg 
c         locfmo(2,1,ibdg) contains the right side fragment for ibdg 
c         locfmo(1,2,ibdg) contains the left side AO for ibdg 
c         locfmo(2,2,ibdg) contains the right side AO for ibdg 
          if(ifmostp.le.2) then
            locfmo(iside+1,1,ibdg)=icurfg
            locfmo(iside+1,2,ibdg)=iloc
          endif
c         if(ifmostp.le.2) locfmo(icurfg,ibdg)=iloc
c         if(ifmostp.eq.1) locfmo(icurfg,ibdg)=iloc
          if(locsav) goto 300
c
c         Rotate HMO LCAO coefficients: find the bond direction 
c
          call vsub(fmoc(1,leftend),1,fmoc(1,rightend),1,bond,1,3)
c
c         now we destroy the first matrix of the point group and
c         store the new rotation matrix (from the standard orientation
c         (along z-axis) to the current bond direction).
c         Then we generate AO transformation matrices.
c
          call vecrot(zaxis,bond,tt)
c         call TRPOSQ(tt,3)
          call trmat
          ibdtyp3=ibdtyp
          if(modQAbas) ibdtyp3=ibdtyp+(iside+1)*maxknd
          NAO=nCBS(ibdtyp)
          nmo=nCAO(ibdtyp3)
          if(modQAbas) call vclr(rotlcao,1,maxcbs*nao) 
c         write(6,*) 'wwwLMOs before rot',ibdtyp3,nmo,nao,maxcbs
c         call prsq(CoreAO(1,1,ibdtyp3),nmo,nao,maxcbs)
          call rotcao(jj,jat,PTR,DTR,FTR,GTR,nao,nmo,CoreAO(1,1,ibdtyp3)
     *              ,maxcbs,rotlcao,maxcbs,nshell,katom,KTYPE,kloc,kmin)
          if(modQAbas) then
           call copyZbl(nao,nmo,l1,iloc,wrk2,rotlcao,maxCBS,ss,l1,nsalc)
          else
c
c         Compute weighted density (based on HO LCAO coefficients)
c
          nao2=(nao*nao+nao)/2
          call vclr(dd,1,nao2)
          ifound=0
          do imo=1,nmo
c           left and right sides have 5 orbitals (1s2s2p) divided.
            if(iside.eq.0.and.iaprjo(imo,ibdtyp).ne.0.or.
     *         iside.ne.0.and.japrjo(imo,ibdtyp).ne.0) then 
              if(some) write(6,*) '  Found i-HMO',imo
              ifound=ifound+1
              iloop=1
              bshift=shiftb(imo,ibdtyp)
              if(.not.saveh) bshift=-1.0D+00
c             bshift=shiftb(imo,ibdtyp)*two
              do i=1,nao
                call daxpy(i,bshift*rotlcao(i,imo),rotlcao(1,imo),1,
     *                     dd(iloop),1)
c               call daxpy(i,bshift*CoreAO(i,imo,ibdtyp),
c    *                     CoreAO(1,imo,ibdtyp),1,dd(iloop),1)
                iloop=iloop+i
c               do j=1,i
c                 iloop=iloop+1
c                 dd(iloop)=dd(iloop)+bshift*CoreAO(i,imo,ibdtyp)*
c    *                                       CoreAO(j,imo,ibdtyp)
c               enddo 
              enddo 
            endif
          enddo 
          norbproj=norbproj+ifound
c         write(6,*) 'wwwhereproj',norbproj,ifound
          if(some) then
            write(6,*) 'HO weighted density, MO '
            call prtrile(dd,nao)
          endif
c
c         compute contribution from this bond projector: S-t * D * S.
c         and add it to the Hamiltonian.
c           
          if(saveh) then
            call TFTRI(wrk1,dd,ss(iloc,1),WRK2,l1,nao,l1)
            if(some) then
              write(6,*) 'Contribution to the Fock matrix from P',ibdg
              call prtrile(wrk1,l1)
            endif
            if(ihnew.eq.0.and.saveh.and..not.parrgo.and..not.modQMbas)
     *         CALL DAread(IDAF,IODA,H,L2,11,0)
            ihnew=ihnew+ifound
            call daxpy(l2,one,wrk1,1,h,1)
          else
c           call vclr(wrk1,1,l2)
c           call tvadd(nao,l1,iloc,iloc,dd,1,wrk1,1)
c           call tritri('c',nao,nao,l1,1,1,dd,1,iloc,iloc,wrk1,1)
c           construct 1-dd*S
c           call madtrap(nao,wrk1,l1,iloc,s,ss)
            call madtrap(nao,dd,l1,iloc,s,ss)
          endif
          endif
        endif
  300 continue
      if(parrgo) then
        if(nxt) CALL DDI_DLBRESET
      endif
      if(locsav) return
      if(parrgo) then 
        if(saveh) then
          call ddi_gsumf(2418,h,l2)
          if(.not.modQMbas) ihnew=1
c         enforce writing in parallel to avoid having to decide if to write
        else
          call ddi_gsumf(2418,ss,l1*l1)
        endif
        call ddi_gsumi(2418,norbproj,1)
      endif
      if(modQAbas) then
        call pushZback(ss,wrk2,wrk1,l1)
      else
        if(saveh) then 
          if(modQMbas) then
c           if ihnew is 0, then no projection operator was constructed.
c           This is sometimes possible (see a few lines below).
c           If this is the case, then the Q matrix is generated as usual.
            call FMOqmt(s,q,ss,wrk1,h,scr,wrk2,L1,L2,L3,ihnew,nqmt)
          else
            if(ihnew.ne.0) then
c             write out the Hamiltonian
c             It is possible (but only if separate molecules are present
c             as parts of the supermolecule) that
c             the Hamiltonian is not changed for a given fragment here.
              CALL DAwrit(IDAF,IODA,h,L2,11,0)
            endif
          endif
        else 
c         add the unit matrix to S 
          do i=1,l1
            ss(i,i)=ss(i,i)+one 
          enddo
c         call TFTRI(wrk1,h,ss,WRK2,l1,l1,l1)
c         call prsq(ss,l1,l1,l1)
          call MRARBR(SS,l1,l1,l1,h,l1,l1,wrk1,l1)
          call dcopy(l1*l1,wrk1,1,h,1)
        endif
      endif
      if(modQAbas) then
c       store the new SALC. 
        call dawrit(IDAF,IODA,ss,L3,44,0)
        nsalc2=(nsalc*nsalc+nsalc)/2
        nsalc3=nsalc*nsalc
        IPA(1)=nsalc
c       These assignments brutally assume C1 symmetry. 
        write(6,*) 'wwwadjusted IPA',nsalc,nsafmo,nqmt
      endif
      if(modQMbas) then
c       store the new Q matrix.
        call dawrit(IDAF,IODA,ss,L3,45,0)
        write(6,*) 'wwwadjusted IPA',nsalc,nsafmo,nqmt
      endif
c
c     restore the rotation matrices for the unit matrix 
c
      call RUNITV(3,3,tt)
      call trmat
      return
      END
c
C*MODULE fmoint  *DECK rotcao
      subroutine rotcao(jj,jat,PTR,DTR,FTR,GTR,l1,nmo,clcao,ldc,rotlcao,
     *                  ldr,nshell,katom,KTYPE,kloc,kmin)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical some,GOPARR,DSKWRK,MASWRK
      dimension ptr(3,3),dtr(6,6),ftr(10,10),gtr(15,15),
     *          clcao(ldc,nmo),rotlcao(ldr,nmo),
     *          katom(*),KTYPE(*),kloc(*),kmin(*) 
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      data dbgfmo/8HDBGFMO  /,dbgme/8HROTCAO  /,debug/8HDEBUG   /
c
c    Rotate LCAO coefficients by using precomputed AO transformation matrices
c    (set up in TRMAT, for a given rotation). Rotate AOs starting from shell jj
c    and only on atom jat. if jat is zero then rotate AOs on all atoms.
c    L1 gives the number of AOs on jat and NMO the number of MOs.
c     parstat: NodeNone 
c
      some=(exetyp.eq.debug.or.exetyp.eq.dbgfmo.or.exetyp.eq.dbgme).and.
     *     maswrk
      if(some) then
        write(iw,*) 'orig LCAO'
        call prsq(clcao,nmo,l1,ldc)
      endif
      loc0=0
      if(jj.ne.0) loc0=kloc(jj)-1
      nao=0
      do ii=jj,nshell
        if(katom(ii).ne.jat.and.jat.ne.0) goto 100 
        mini=kmin(ii)
        loci=kloc(ii)-loc0
        LIT = KTYPE(II) 
c       S and the S part of L: spherically symmetric
        if(mini.eq.1) then
          call dcopy(nmo,clcao(loci,1),ldc,rotlcao(loci,1),ldr)
          nao=nao+1
c         shift to P for L shells
          if(lit.eq.2) loci=loci+1
        endif
        if(lit.eq.2) then
          call MRARBR(ptr, 3, 3, 3,clcao(loci,1),ldc,nmo,rotlcao(loci,1)
     *               ,ldr)
          nao=nao+3
        else if(lit.eq.3) then 
          call MRARBR(dtr, 6, 6, 6,clcao(loci,1),ldc,nmo,rotlcao(loci,1)
     *               ,ldr)
          nao=nao+6
        else if(lit.eq.4) then
          call MRARBR(ftr,10,10,10,clcao(loci,1),ldc,nmo,rotlcao(loci,1)
     *               ,ldr)
          nao=nao+10
        else if(lit.eq.5) then
          call MRARBR(gtr,15,15,15,clcao(loci,1),ldc,nmo,rotlcao(loci,1)
     *               ,ldr)
          nao=nao+15
        endif
      enddo
  100 continue
      if(nao.ne.l1) then
c       not all shells found in the basis set?! 
        if(maswrk) write(iw,9000) nao,l1
        call abrt
      endif
      if(some) then
        write(iw,*) 'rot LCAO'
        call prsq(rotlcao,nmo,l1,ldr)
      endif
      return
 9000 format(/1x,'Basis set size mismatch in ROTCAO:',2I5,
     *       /1x,'This is probably caused my a mismatch in your',
     *           'basis set and $FMOLMO',
     *     /1X,'(e.g. you use 6-31G* but defined STO-3G in $FMOLMO).',/)

      END
C*MODULE fmoint  *DECK diminid
      SUBROUTINE diminid(itype,ifg,jfg,l1i,l1j,nai,naj,l1d,da,db,wrk,
     *                   wrk1,iodfmo,irec0,iabdfg,jabdfg,indat,iaglob,
     *                   locfmo,mapi,mapj,orbxch,enexch,iodexch,jodexch,
     *                   some)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXSH=5000,MXGTOT=20000,ZERO=0.0D+00,two=2.0D+00)
      logical isini,isinj,ok,ghostini,ghostinj,orbxch,enexch,iodexch,
     *        jodexch,some,dbg
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),ch(mxgtot),ci(mxgtot),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmorun/ espscf,IDAFMO,icurlay,icurfg,jcurfg,kcurfg,ncursh,
     *                ngau,icurpop,ifmostp,moncor,nat1e,needr,modrst,
     *                norbproj
      dimension da(*),db(*),wrk(*),wrk1(*),iodfmo(*),iabdfg(*),jabdfg(*)
     *         ,indat(*),iaglob(*),locfmo(2,2,*),mapi(*),mapj(*)
c
c     Prepare initial density for a dimer (including BSSE case)
c     itype=0 add two monomer densities
c     itype=1 use ifg monomer density and zero for monomer jfg
c     itype=2 use jfg monomer density and zero for monomer ifg
c     itype=-1 light version: draw maps and quit (no density made) 
c     If modorb.ne.0 then instead of density first construct the orbitals,
c     save them and then build density.
c     parstat: GroupNone 
c
c     DB/DA are used as temporary storage
c
      dbg=.false.
      l2i=(l1i*l1i+l1i)/2
      l2j=(l1j*l1j+l1j)/2
      l2d=(l1d*l1d+l1d)/2
      l3d=l1d*l1d
      l3i=l1i*l1i
      l3j=l1j*l1j
      idaoff=l2d
      if(orbxch) idaoff=idaoff+l3d
      if(orbxch) then 
        m2i=l3i
        m2j=l3j
      else
        m2i=l2i
        m2j=l2j
      endif
      if(enexch) then
        m2i=m2i+l1i
        m2j=m2j+l1j
      endif
      if(orbxch.and.itype.ge.0) call vclr(wrk1,1,l1d)
      if(itype.gt.0) then
        m2d=l2d
        if(orbxch) m2d=m2d+l3d
        if(enexch) m2d=m2d+l1d
        call vclr(da,1,m2d)
        goto 200
      endif
      call viclr(mapi,1,l1d)
      call viclr(mapj,1,l1d)
      loopi=0
      loopj=0
      loopij=0
      iextra=0
      iextraa=0
      jextra=0
      jextraa=0
      iatprev=-1
      do i=1,nshell
        iat=katom(i)
        iatg=iaglob(iat)
        mini=kmin(i)
        maxi=kmax(i)
        ijfg=indat(iatg)
        ok=.false.
        loopij0=loopij
        if(iatprev.ne.iat) then
          iextraa=iextraa+iextra
          iextra=0
          jextraa=jextraa+jextra
          jextra=0
          iatprev=iat
c         This is needed to reset index for ghost atoms (the index 
c         runs continuosly within a given atom, and it is reset when 
c         the atom iat changes. jextra runs for one atom; jextraa 
c         accumulates the number of added AOs due to ghost atoms; 
c         this variable is used only for double checking.
        endif
c       ibdfg=0 corresponds to handling atoms in fragments themselves
        do 100 ibdfg=0,nbdfg
          if(ibdfg.ne.0) then
            iatb=iabdfg(ibdfg)
            jatb=jabdfg(ibdfg)
c           the job aborts because of the present limitation of the code
c           in this subroutine; it should be possible to fix if needed.
            if(iatb.ge.0.or.jatb.le.0) call abrt
            iatb=abs(iatb)
            jatb=abs(jatb)
            jfrgb=indat(jatb)
            ghostini=iatb.eq.iatg.and.jfrgb.eq.ifg
            ghostinj=iatb.eq.iatg.and.jfrgb.eq.jfg
            isini=ghostini
            isinj=ghostinj
            iloc=0
            jloc=0
            do iside=1,2
              locfg=locfmo(iside,1,ibdfg)
              locao=locfmo(iside,2,ibdfg)
              if(locfg.eq.ifg) iloc=locao
              if(locfg.eq.jfg) jloc=locao
c             It is not possible that both sides are in the same fragment
c             (with sane input).
            enddo 
c           iloc=locfmo(ifg,ibdfg) 
c           jloc=locfmo(jfg,ibdfg) 
          else
            ghostini=.false.
            ghostinj=.false.
            isini=ijfg.eq.ifg
            isinj=ijfg.eq.jfg
            iloc=0
            jloc=0
          endif
          if(.not.(isini.or.isinj)) goto 100
          loopij=loopij0
          do ii=mini,maxi
            loopij=loopij+1
            if(isini) then
              if(iloc.ne.0.and.ghostini) then
                mapi(loopij)=iloc+iextra
                iextra=iextra+1
              else
                loopi=loopi+1
                mapi(loopij)=loopi
              endif
              ok=.true.
            endif
            if(isinj) then
cl            if(iloc.ne.0.and.jloc.ne.0.and.ghostinj) then
              if(jloc.ne.0.and.ghostinj) then
                mapj(loopij)=jloc+jextra
                jextra=jextra+1
              else
                loopj=loopj+1
                mapj(loopij)=loopj
              endif
              ok=.true.
            endif
          enddo
c         if(ok) goto 110
  100   continue
        if(.not.ok) then
          write(6,*) 'Atom not found in diminid',iat,iaglob(iat)
          call abrt
        endif
      enddo
c     add up whatever left from added orbitals.
      iextraa=iextraa+iextra
      jextraa=jextraa+jextra
      if(loopi+iextraa.ne.l1i.or.loopj+jextraa.ne.l1j.or.loopij.ne.l1d)
     *  then
        write(6,*) 'Confusion in diminid:',loopi,iextraa,l1i,loopj,
     *                                       jextraa,l1j,loopij,l1d
        call abrt
      endif
  200 continue
      if(some) then
        write(6,*) 'Final mapi:',(mapi(i),i=1,l1d)
        write(6,*) 'Final mapj:',(mapj(i),i=1,l1d)
      endif
      if(itype.lt.0) return
      call mapcheck(mapi,wrk,l1i,l1d,imap)
      call mapcheck(mapj,wrk,l1j,l1d,jmap)
      if(imap.ne.0.or.jmap.ne.0) then
         write(6,*) 'maps are false!',imap,jmap
         call abrt
      endif
      if(itype.ne.2) then
        m2i=l2i
        if(orbxch) m2i=l3i
        ioff=0
        if(orbxch) ioff=l2i
        ioffe=ioff+m2i
        if(enexch) m2i=m2i+l1i
        if(iodexch) then
          CALL rareads(IDAFMO,iodfmo,wrk,l2i+m2i,ifg+irec0,0)
        else
          CALL rareads(IDAFMO,iodfmo,wrk(ioff+1),m2i,ifg+irec0,0)
          if(orbxch) call DMTX2(wrk,wrk(ioff+1),nai,l1i,l1i) 
        endif
      endif 
      if(dbg) then
        write(6,*) 'wwwDensity I',ifg
        call prtril(wrk,l1i) 
        if(orbxch) then 
          write(6,*) 'wwworbital I',nai,ioff,m2i
          call prsq(wrk(ioff+1),l1i,l1i,l1i)
        endif 
        if(enexch) write(6,*) 'wwweni',(wrk(ioffe+ii),ii=1,l1i)
      endif
      if(itype.ne.1) then
        m2j=l2j
        if(orbxch) m2j=l3j
        joff=0
        if(orbxch) joff=l2j
        joffe=joff+m2j
        if(enexch) m2j=m2j+l1j
        if(jodexch) then
          CALL rareads(IDAFMO,iodfmo,db,l2j+m2j,jfg+irec0,0)
        else
          CALL rareads(IDAFMO,iodfmo,db(joff+1),m2j,jfg+irec0,0)
          if(orbxch) call DMTX2(db,db(joff+1),naj,l1j,l1j)
        endif
      endif
      if(dbg) then
        write(6,*) 'wwwDensity J',jfg
        call prtril(db,l1j)
        if(orbxch) then
          write(6,*) 'wwworbital J',naj
          call prsq(db(joff+1),l1j,l1j,l1j)
        endif
        if(enexch) write(6,*) 'wwwenj',(db(joffe+ii),ii=1,l1j)
      endif
c     after all this bizarreness, wrk and db will hold:
c     D+V+E (in this order), any of V and E can be omitted if corresponding
c     options (orbxch and enexch) are not set. wrk is frag I and db frag J.
c     may be safer to fill high energies
      if(enexch) call vclr(da(idaoff+1),1,l1d)
      if(enexch.and..not.orbxch)
     *  write(6,*) 'Warning: orbital energies are unreliable.'
c     tol should be set approximately 2-3 times smaller than the smalles 
c     shiftb
c     tol=1.0D+05
      loop=0
      loop1=l2d
c     the loop order i,j is important for orbxch=.true. (see loop)
      do i=1,l1d
        jend=i
        if(orbxch) jend=l1d
c       due to the loop structure for orbxch=.t. runs i corresponds to
c       MO and j to AO and not otherwise (contrary to normal case).
c
c       imon and jmon keep the number of LCAO coefficients taken from
c       ifg and jfg for the MO i. It helps deciding which orbital energy
c       to use. This information is passed through array da at the location
c       where dimer energies are stored. Pristine value (ZERO) means
c       use I energies, otherwise use J. 
c       Trouble happens if we only have D and E (not V). Then no reliable
c       way is implemented for orbital energies. This combination is
c       at present of no use, though. 
c       
c       imon=0
c       jmon=0
        do j=1,jend
          ii=mapi(i)
          ij=mapi(j)
          ji=mapj(i)
          jj=mapj(j)
          if(i.ge.j) loop=loop+1
          loop1=loop1+1
          ii1=max(ii,ij)
          ij1=min(ii,ij)
          ji1=max(ji,jj)
          jj1=min(ji,jj)
          isini=ii.ne.0.and.ij.ne.0.and.itype.ne.2
          isinj=ji.ne.0.and.jj.ne.0.and.itype.ne.1 
c         if(ii.ne.0.and.ij.ne.0.and.itype.ne.2) then
          if(isini) then
c    *       (isinj .or. ji.eq.0.and.jj.eq.0.and.itype.ne.2)) then
            if(orbxch) then
c             this enforces using MO coefficients only for I-occupied
c             orbitals or if J-orbitals are not defined.
c             the last condition ignores orbitals with lunatic energies
c             Can it generate zero orbitals?? (if neither I nor J is used)
c             if(ii.le.nai) then
              if((ii.le.nai.or.ji.eq.0.or.ji.gt.naj)) then
c    *           .and.wrk(ioffe+ii).lt.tol) then
                da(loop1)=wrk(ioff+(ii-1)*l1i+ij)
c               imon=imon+1
              else
                da(loop1)=zero
              endif
            endif
c           Fill in density only for lower triangle of indices
            if(i.ge.j) da(loop)=wrk((ii1*ii1-ii1)/2+ij1)
          else
            if(orbxch) da(loop1)=zero
            if(i.ge.j) da(loop)=zero
          endif
c         if(ji.ne.0.and.jj.ne.0.and.itype.ne.1) then
          if(isinj) then
c    *       (isini .or. ii.eq.0.and.ij.eq.0.and.itype.ne.1)) then
            if(orbxch) then
              if((ji.le.naj.or.ii.eq.0.or.ii.gt.nai)) then
                da(loop1)=da(loop1)+db(joff+(ji-1)*l1j+jj)
c               jmon=jmon+1
              endif
            endif
            if(i.ge.j) 
     *        da(loop)=da(loop)+db((ji1*ji1-ji1)/2+jj1)
          endif
        enddo
c       if(orbxch.and.jmon.gt.imon) da(idaoff+i)=two 
c       if(orbxch) write(6,*) i,'-th MO, ',imon,jmon
c       anything nonzero will do.
      enddo
c     fill in occupation numbers and orbital energies
      if(orbxch.or.enexch) then
        do i=1,l1d
          ii=mapi(i)
          ji=mapj(i)
c         ei=da(idaoff+i)
          if(ii.ne.0.and.itype.ne.2) then
            if(orbxch.and.ii.le.nai) wrk1(i)=two
            if(enexch) then
c             if(orbxch) then
c               "reliable" way based on LCAO copy counter
c               if(ei.eq.zero) da(idaoff+i)=wrk(ioffe+ii)
c             else
c               unreliable, a better way is to use nai/naj to aid.
                da(idaoff+i)=wrk(ioffe+ii)
c             endif
            endif
          endif
          if(ji.ne.0.and.itype.ne.1) then
c           it seems possible that wrk1(i) becomes equal to 4. But in
c           physically reasonable systems it should not happen.
            if(orbxch.and.ji.le.naj) wrk1(i)=wrk1(i)+two
            if(enexch) then
c             if(orbxch) then
c               "reliable" way based on LCAO copy counter
c               if(ei.ne.zero) da(idaoff+i)=db(joffe+ji)
c             else
                if(da(idaoff+i).eq.zero) then
                  da(idaoff+i)=db(joffe+ji)
                else
                  da(idaoff+i)=min(da(idaoff+i),db(joffe+ji))
                endif
c               da(idaoff+i)=da(idaoff+i)+db(joffe+ji)
c               should we instead do averaging?
c             endif
            endif
          endif
        enddo
      endif
      if(orbxch) then
         if((iodexch.or.jodexch).and.itype.eq.0) then
           idum=0
           call fmoord(da(l2d+1),da(idaoff+1),iodexch,jodexch,.false.,
     *                 mapi,mapj,idum,enexch,l1d,nai,naj,0,wrk,wrk1)
         else
           if(enexch) call ORDERV(da(l2d+1),da(idaoff+1),wrk,wrk1,l1d,
     *                            l1d,l1d)
         endif
         CALL dawrit(IDAF,IODA,da(l2d+1),L3d,15,0)
         if(dbg) then
           write(6,*) 'wwworbital4-final'
           call prsq(da(l2d+1),l1d,l1d,l1d)
         endif
      endif
      if(enexch) CALL dawrit(IDAF,IODA,da(idaoff+1),L1d,17,0)
      CALL dawrit(IDAF,IODA,da,L2d,16,0)
      if(dbg) then
        if(enexch) write(6,*) 'wwwen',(da(idaoff+ii),ii=1,l1d)
        write(6,*) 'wwwDensity4-final'
        call prtril(da,l1d)
      endif
      IF(some) THEN
         WRITE(IP,*) 'Orbitals cooked up by diminid for',ifg,jfg
         WRITE(IP,*) '$VEC'
         CALL PUSQL(da(l2d+1),l1d,l1d,l1d)
         WRITE(IP,*) '$END'
      END IF
      RETURN
      END
C*MODULE fmoint  *DECK mapcheck
      subroutine mapcheck(map,iwrk,n,m,ires)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      dimension map(m),iwrk(n)
c
c     check if array map contains all integers from 1 to n, without 
c     redundancies.
c     Redundant 0 elements are allowed and skipped when checking.
c     return 0 if each integer occurs only once and the offending integer 
c     otherwise 
c              -1 if map contains numbers larger than n or negative integers
c
      ires=-1
      call viclr(iwrk,1,n)
      do i=1,m
        mapi=map(i)
        if(mapi.gt.n.or.mapi.lt.0) return
        if(mapi.ne.0) iwrk(mapi)=iwrk(mapi)+1
      enddo
      do i=1,n
        if(iwrk(i).ne.1) then
          ires=i
          return
        endif
      enddo
      ires=0
      RETURN
      END
C*MODULE fmoint  *DECK mod2ei
      subroutine mod2ei 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000)
      logical PK,PACK2E,PANDK,BLOCK
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      common /INT2IC/ NINTIC,iNINTIC,nxxic,lbufpic,lixic,labsix,nintix
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,inttyp,igrdtyp
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /PCKLAB/ LABSIZ
      COMMON /PKFIL / PK,PANDK,BLOCK
c
c     if we constructed the basis set we should readjust LABSIZ
c     and do memory reallocation. This relies on coherence with
c     int2eic in INT2A.SRC. Caveat emptor!
c     LABSIZ may be changed frantically during FMO runs!
c     Below the only allowed values for LABSIZ is 1 and 2.
c
      MAXAO=255
      IF(PK) MAXAO=361
      LABSIZ = 1
c     IF(NUM.GT.MAXAO.and.nwdvar.eq.2) LABSIZ = 2
      IF(NUM.GT.MAXAO) LABSIZ = 2
c     Check if the user set LABSIZ incorrectly.
      if(labsix.lt.LABSIZ) then
         write(iw,9000) labsiz,labsix
         call abrt 
      endif
      if(NINTIC.eq.0) return
c     labsiz2=2/labsiz
      NINTIC=NINTIx
c     recompute NINTIC if initially labsiz was set to 2.
c     There is nothing to be done if labsiz=labsix.
c     That is, the only adjustment is needed from LABSIZ=2 to LABSIZ=1.
      if(labsiz.lt.labsix) then
c       The formula is obtained as follows:
c       (2 in the RHS is 2/labsiz' and in the LHS 2/labsiz is omitted
c       as the only valid combination is labsiz'=1 and labsiz=2). 
c       real  int    real   integer
c       N+M + N+M = N'+M + (N'+M)/2
c       N=NINTIx, M=nINTMX, N=NINTIC, and LHS&RHS are equal to NEED(below)
c       LHS: NEED with original LABSIZ (LABSIX)
c       RHS: NEED with reduced labsiz. 
c       nINTMX+NINTIx should be even for the formula to work or else
c       substract one(?)
c       Round off should be from below!
        NINTIC=(nINTMX+4*NINTIx)/3
c       write(6,*) 'Warning: increasing NINTIC to ',NINTIC,NINTIx
      endif 
      iNINTIC=NINTIC
      if(labsiz.eq.2.and.nwdvar.eq.2) iNINTIC=NINTIC*2
      if(labsiz.eq.1.and.nwdvar.eq.1) iNINTIC=NINTIC/2
c     Reallocate the memory!
c     Since even presumably unloaded guns sometimes kill, do this check.
      if(lbufpic.eq.0) call abrt
c     The following must be used in INT2A.SRC: 
c     lbufpic= LOADFM + 1
c     lixic=   lbufpic+ nINTMX+NINTIC
c     LAST=    lixic  + (nINTMX+NINTIC-1)/labsiz2+1
c     NEED=    LAST - LOADFM - 1
c     LABSIZ2=2/LABSIZ
c     newneed=nINTMX+NINTIC+(nINTMX+NINTIC-1)/labsiz2+1
c
      lixic=   lbufpic+ nINTMX+NINTIC
c     write(6,9010) LABSIZ,NINTIC,ININTIC,newneed
c
      return
 9000 format(/1x,'LABSIZ needs to be',I2,
     *           ' whereas it is set in $INTGRL to',I2,/)
c9010 format(/1x,'Int dimensions: LABSIZ=',I1,', NINTIC=',I9,
c    *           ' ININTIC=',I9,',NEED=',I10)
      end
C*MODULE fmoint  *DECK triminid
      SUBROUTINE triminid(ifg,jfg,kfg,l1i,l1j,l1k,nai,naj,nak,l1t,da,db,
     *                    dc,wrk,wrk1,iodfmo,irec0,iabdfg,jabdfg,indat,
     *                    iaglob,locfmo,mapi,mapj,mapk,orbxch,enexch,
     *                    iodexch,jodexch,kodexch,some)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXSH=5000,MXGTOT=20000,ZERO=0.0D+00,two=2.0D+00)
      logical isini,isinj,isink,ok,ghostini,ghostinj,ghostink,orbxch,
     *        enexch,iodexch,jodexch,kodexch,some
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),ch(mxgtot),ci(mxgtot),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmorun/ espscf,IDAFMO,icurlay,icurfg,jcurfg,kcurfg,ncursh,
     *                ngau,icurpop,ifmostp,moncor,nat1e,needr,modrst,
     *                norbproj
      dimension da(*),db(*),dc(*),wrk(*),wrk1(*),iodfmo(*),iabdfg(*),
     *          jabdfg(*),indat(*),iaglob(*),locfmo(2,2,*),mapi(*),
     *          mapj(*),mapk(*)
c
c     Prepare initial density for a trimer
c     itype=0 add two monomer densities
c     If modorb.ne.0 then instead of density first construct the orbitals,
c     save them and then build density.
c     Note that at present only the density construction works.
c     parstat: GroupNone 
c
c     DB/DA are used as temporary storage
c
      l2i=(l1i*l1i+l1i)/2
      l2j=(l1j*l1j+l1j)/2
      l2k=(l1k*l1k+l1k)/2
      l2t=(l1t*l1t+l1t)/2
      l3t=l1t*l1t
      l3i=l1i*l1i
      l3j=l1j*l1j
      l3k=l1k*l1k
      idaoff=l2t
      if(orbxch) idaoff=idaoff+l3t
      if(orbxch) then 
        m2i=l3i
        m2j=l3j
        m2k=l3k
      else
        m2i=l2i
        m2j=l2j
        m2k=l2k
      endif
      if(enexch) then
        m2i=m2i+l1i
        m2j=m2j+l1j
        m2k=m2k+l1k
      endif
      if(orbxch) call vclr(wrk1,1,l1t)
      call viclr(mapi,1,l1t)
      call viclr(mapj,1,l1t)
      call viclr(mapk,1,l1t)
      loopi=0
      loopj=0
      loopk=0
      loopij=0
      iextra=0
      iextraa=0
      jextra=0
      jextraa=0
      kextra=0
      kextraa=0
      iatprev=-1
      do i=1,nshell
        iat=katom(i)
        iatg=iaglob(iat)
        mini=kmin(i)
        maxi=kmax(i)
        ijfg=indat(iatg)
        ok=.false.
        loopij0=loopij
        if(iatprev.ne.iat) then
          iextraa=iextraa+iextra
          iextra=0
          jextraa=jextraa+jextra
          jextra=0
          kextraa=kextraa+kextra
          kextra=0
          iatprev=iat
c         This is needed to reset index for ghost atoms (the index 
c         runs continuosly within a given atom, and it is reset when 
c         the atom iat changes. jextra runs for one atom; jextraa 
c         accumulates the number of added AOs due to ghost atoms; 
c         this variable is used only for double checking.
        endif
c       ibdfg=0 corresponds to handling atoms in fragments themselves
        do 100 ibdfg=0,nbdfg
          if(ibdfg.ne.0) then
            iatb=iabdfg(ibdfg)
            jatb=jabdfg(ibdfg)
c           the job aborts because of the present limitation of the code
c           in this subroutine; it should be possible to fix if needed.
            if(iatb.ge.0.or.jatb.le.0) call abrt
            iatb=abs(iatb)
            jatb=abs(jatb)
            jfrgb=indat(jatb)
            ghostini=iatb.eq.iatg.and.jfrgb.eq.ifg
            ghostinj=iatb.eq.iatg.and.jfrgb.eq.jfg
            ghostink=iatb.eq.iatg.and.jfrgb.eq.kfg
            isini=ghostini
            isinj=ghostinj
            isink=ghostink
            iloc=0
            jloc=0
            kkloc=0
            do iside=1,2
              locfg=locfmo(iside,1,ibdfg)
              locao=locfmo(iside,2,ibdfg)
              if(locfg.eq.ifg) iloc=locao
              if(locfg.eq.jfg) jloc=locao
              if(locfg.eq.kfg) kkloc=locao
c             It is not possible that both sides are in the same fragment
c             (with sane input).
            enddo 
          else
            ghostini=.false.
            ghostinj=.false.
            ghostink=.false.
            isini=ijfg.eq.ifg
            isinj=ijfg.eq.jfg
            isink=ijfg.eq.kfg
            iloc=0
            jloc=0
            kkloc=0
          endif
          if(.not.(isini.or.isinj.or.isink)) goto 100
          loopij=loopij0
          do ii=mini,maxi
            loopij=loopij+1
            if(isini) then
              if(iloc.ne.0.and.ghostini) then
                mapi(loopij)=iloc+iextra
                iextra=iextra+1
              else
                loopi=loopi+1
                mapi(loopij)=loopi
              endif
              ok=.true.
            endif
            if(isinj) then
              if(jloc.ne.0.and.ghostinj) then
                mapj(loopij)=jloc+jextra
                jextra=jextra+1
              else
                loopj=loopj+1
                mapj(loopij)=loopj
              endif
              ok=.true.
            endif
            if(isink) then
              if(kkloc.ne.0.and.ghostink) then
                mapk(loopij)=kkloc+kextra
                kextra=kextra+1
              else
                loopk=loopk+1
                mapk(loopij)=loopk
              endif
              ok=.true.
            endif
          enddo
c         if(ok) goto 110
  100   continue
        if(.not.ok) then
          write(6,*) 'Atom not found in diminid',iat,iaglob(iat)
          call abrt
        endif
      enddo
c     add up whatever left from added orbitals.
      iextraa=iextraa+iextra
      jextraa=jextraa+jextra
      kextraa=kextraa+kextra
      if(loopi+iextraa.ne.l1i.or.loopj+jextraa.ne.l1j.or.
     *   loopk+kextraa.ne.l1k.or.loopij.ne.l1t) then
        write(6,*) 'Confusion in triminid:',loopi,iextraa,l1i,loopj,
     *             jextraa,l1j,loopk,kextraa,l1k,loopij,l1t
        call abrt
      endif
c 200 continue
      if(some) then
        write(6,*) 'Final mapi:',(mapi(i),i=1,l1t)
        write(6,*) 'Final mapj:',(mapj(i),i=1,l1t)
        write(6,*) 'Final mapk:',(mapk(i),i=1,l1t)
      endif
      call mapcheck(mapi,wrk,l1i,l1t,imap)
      call mapcheck(mapj,wrk,l1j,l1t,jmap)
      call mapcheck(mapk,wrk,l1k,l1t,kmap)
      if(imap.ne.0.or.jmap.ne.0.or.kmap.ne.0) then
         write(6,*) 'maps are false!',imap,jmap,kmap
         call abrt
      endif
c     if(itype.ne.2) then
        m2i=l2i
        if(orbxch) m2i=l3i
        ioff=0
        if(orbxch) ioff=l2i
        ioffe=ioff+m2i
        if(enexch) m2i=m2i+l1i
        if(iodexch) then
          CALL rareads(IDAFMO,iodfmo,wrk,l2i+m2i,ifg+irec0,0)
        else
          CALL rareads(IDAFMO,iodfmo,wrk(ioff+1),m2i,ifg+irec0,0)
          if(orbxch) call DMTX2(wrk,wrk(ioff+1),nai,l1i,l1i) 
        endif
        m2j=l2j
        if(orbxch) m2j=l3j
        joff=0
        if(orbxch) joff=l2j
        joffe=joff+m2j
        if(enexch) m2j=m2j+l1j
        if(jodexch) then
          CALL rareads(IDAFMO,iodfmo,db,l2j+m2j,jfg+irec0,0)
        else
          CALL rareads(IDAFMO,iodfmo,db(joff+1),m2j,jfg+irec0,0)
          if(orbxch) call DMTX2(db,db(joff+1),naj,l1j,l1j)
        endif
        m2k=l2k
        if(orbxch) m2k=l3k
        koff=0
        if(orbxch) koff=l2k
        koffe=koff+m2k
        if(enexch) m2k=m2k+l1k
        if(kodexch) then
          CALL rareads(IDAFMO,iodfmo,dc,l2k+m2k,kfg+irec0,0)
        else
          CALL rareads(IDAFMO,iodfmo,dc(koff+1),m2k,kfg+irec0,0)
          if(orbxch) call DMTX2(dc,dc(koff+1),nak,l1k,l1k)
        endif
c     after all this bizarreness, wrk and db will hold:
c     D+V+E (in this order), any of V and E can be omitted if corresponding
c     options (orbxch and enexch) are not set. wrk is frag I and db frag J.
c     may be safer to fill high energies
      if(enexch) call vclr(da(idaoff+1),1,l1t)
      if(enexch.and..not.orbxch)
     *  write(6,*) 'Warning: orbital energies are unreliable.'
      loop=0
      loop1=l2t
c     the loop order i,j is important for orbxch=.true. (see loop)
      do i=1,l1t
        jend=i
        if(orbxch) jend=l1t
c       due to the loop structure for orbxch=.t. runs i corresponds to
c       MO and j to AO and not otherwise (contrary to the convention).
c
        do j=1,jend
          ii=mapi(i)
          ij=mapi(j)
          ji=mapj(i)
          jj=mapj(j)
          ki=mapk(i)
          kj=mapk(j)
          if(i.ge.j) loop=loop+1
          loop1=loop1+1
          ii1=max(ii,ij)
          ij1=min(ii,ij)
          ji1=max(ji,jj)
          jj1=min(ji,jj)
          ki1=max(ki,kj)
          kj1=min(ki,kj)
c         isini=ii.ne.0.and.ij.ne.0.and.itype.ne.2
c         isinj=ji.ne.0.and.jj.ne.0.and.itype.ne.1 
          isini=ii.ne.0.and.ij.ne.0
          isinj=ji.ne.0.and.jj.ne.0
          isink=ki.ne.0.and.kj.ne.0
c         if(ii.ne.0.and.ij.ne.0.and.itype.ne.2) then
          if(isini) then
c    *       (isinj .or. ji.eq.0.and.jj.eq.0.and.itype.ne.2)) then
            if(orbxch) then
c             this enforces using MO coefficients only for I-occupied
c             orbitals or if J-orbitals are not defined.
c             the last condition ignores orbitals with lunatic energies
c             Can it generate zero orbitals?? (if neither I nor J is used)
              if(ii.le.nai.or.
     *         ((ji.eq.0.or.ji.gt.naj).and.(ki.eq.0.or.ki.gt.nak))) then
                da(loop1)=wrk(ioff+(ii-1)*l1i+ij)
              else
                da(loop1)=zero
              endif
            endif
c           Fill in density only for lower triangle of indices
            if(i.ge.j) da(loop)=wrk((ii1*ii1-ii1)/2+ij1)
          else
            if(orbxch) da(loop1)=zero
            if(i.ge.j) da(loop)=zero
          endif
c         if(ji.ne.0.and.jj.ne.0.and.itype.ne.1) then
          if(isinj) then
c    *       (isini .or. ii.eq.0.and.ij.eq.0.and.itype.ne.1)) then
            if(orbxch) then
              if(ji.le.naj.or.
     *         ((ii.eq.0.or.ii.gt.nai).and.(ki.eq.0.or.ki.gt.nak))) then
                da(loop1)=da(loop1)+db(joff+(ji-1)*l1j+jj)
              endif
            endif
            if(i.ge.j) 
     *        da(loop)=da(loop)+db((ji1*ji1-ji1)/2+jj1)
          endif
          if(isink) then
c    *       (isini .or. ii.eq.0.and.ij.eq.0.and.itype.ne.1)) then
            if(orbxch) then
              if(ki.le.nak.or.
     *         ((ii.eq.0.or.ii.gt.nai).and.(ji.eq.0.or.ji.gt.naj))) then
                da(loop1)=da(loop1)+dc(koff+(ki-1)*l1k+kj)
              endif
            endif
            if(i.ge.j)
     *        da(loop)=da(loop)+dc((ki1*ki1-ki1)/2+kj1)
          endif
        enddo
      enddo
c     fill in occupation numbers and orbital energies
      if(orbxch.or.enexch) then
        do i=1,l1t
          ii=mapi(i)
          ji=mapj(i)
          ki=mapk(i)
c         if(ii.ne.0.and.itype.ne.2) then
          if(ii.ne.0) then
            if(orbxch.and.ii.le.nai) wrk1(i)=two
            if(enexch) then
c             unreliable, a better way is to use nai/naj to aid.
              da(idaoff+i)=wrk(ioffe+ii)
            endif
          endif
c         if(ji.ne.0.and.itype.ne.1) then
          if(ji.ne.0) then
            if(orbxch.and.ji.le.naj) wrk1(i)=wrk1(i)+two
            if(enexch) then
              if(da(idaoff+i).eq.zero) then
                da(idaoff+i)=db(joffe+ji)
              else
                da(idaoff+i)=min(da(idaoff+i),db(joffe+ji))
              endif
c             should we instead do averaging?
            endif
          endif
          if(ki.ne.0) then
            if(orbxch.and.ki.le.nak) wrk1(i)=wrk1(i)+two
            if(enexch) then
              if(da(idaoff+i).eq.zero) then
                da(idaoff+i)=dc(koffe+ki)
              else
                da(idaoff+i)=min(da(idaoff+i),dc(koffe+ki))
              endif
            endif
          endif
        enddo
      endif
      if(orbxch) then
c        if((iodexch.or.jodexch.or.kodexch).and.itype.eq.0) then
         if(iodexch.or.jodexch.or.kodexch) then
           call fmoord(da(l2t+1),da(idaoff+1),iodexch,jodexch,kodexch,
     *                 mapi,mapj,mapk,enexch,l1t,nai,naj,nak,wrk,wrk1)
         else
           if(enexch) call ORDERV(da(l2t+1),da(idaoff+1),wrk,wrk1,l1t,
     *                            l1t,l1t)
         endif
         CALL dawrit(IDAF,IODA,da(l2t+1),L3t,15,0)
      endif
      if(enexch) CALL dawrit(IDAF,IODA,da(idaoff+1),L1t,17,0)
      CALL dawrit(IDAF,IODA,da,L2t,16,0)
      IF(some) THEN
         WRITE(IP,*) 'Orbitals cooked up by triminid for',ifg,jfg,kfg
         WRITE(IP,*) '$VEC'
         CALL PUSQL(da(l2t+1),l1t,l1t,l1t)
         WRITE(IP,*) '$END'
      END IF
      RETURN
      END
C*MODULE GRD2A   *DECK ESD2der
      SUBROUTINE ESD2der(DI,DJ,l1i)
C
C       THE DRIVER FOR THE TWO ELECTRON GRADIENT contributions to separated
c       dimer energies.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL SKIPI,SKIPJ,SKIPK,SKIPL
      LOGICAL UHFTYP,PACK2E,POPLE,HONDO
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT
      LOGICAL SOME,OUT,DBG
C
      PARAMETER (MXGTOT=20000, MXSH=5000, MXATM=2000, MXAO=8192,
     *           MXFRG=1050, MXFGPT=12000)
C
      COMMON /DERSKP/ IIAT,JJAT,KKAT,LLAT,SKIPI,SKIPJ,SKIPK,SKIPL
      COMMON /DLT   / LAT,LBT,LCT,LDT
      COMMON /DSHLNO/ LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJD,KLD,IJ,KL
      COMMON /DSHLT / RTOL,DTOL,VTOL1,VTOL2,VTOLS,OUT,DBG
      COMMON /FMCOM / X(1)
c     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 /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,inttyp,igrdtyp
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      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 /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SHLBAS/ MAXTYP,MAXNUM
      COMMON /TMVALS/ TI,TX,TIM
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      common /fmorun/ espscf,IDAFMO,icurlay,icurfg,jcurfg,kcurfg,ncursh,
     *                ngau,icurpop,ifmostp,moncor,nat1e,needr,modrst,
     *                norbproj
      dimension di(*),dj(*)
C
      COMMON/DERMEM/IWFN,IXCH,INIJG,IGINT,IFINT,ISINT,IIJKLG,
     1 IDAB,ICHRG,IXY,IXZ,IYZ,IX,IY,IZ,ISJ,ISK,ISL,IGIJKL,IGNKL,IGNM,
     2 IDIJ,IDKL,IB00,IB01,IB10,IC00,ID00,IF00,
     3 IAAI,IAAJ,IBBK,IBBL,IFI,IFJ,IFK,IFL,
     4 ISII,ISJJ,ISKK,ISLL,ISIJ,ISIK,ISIL,ISJK,ISJL,ISKL,
     5 IDIJSI,IDIJSJ,IDKLSK,IDKLSL,IABV,ICV,IRW
      COMMON/INDD80/IMAX,JMAX,KKKMAX,LMAX
C
      PARAMETER (RLN10=2.30258D+00)
      PARAMETER (TEN=10.0D+00, ONE=1.0D+00)
      PARAMETER (TENM9=1.0D-09, TENM11=1.0D-11)
      PARAMETER (TENM20=1.0D-20, PT5=0.5D+00, TENM12=1.0D-12)
C
      DIMENSION LENSHL(5)
      DATA LENSHL/1,4,10,20,35/
      DATA CHECK,GRD2,DEBUG/8HCHECK   ,8HGRD2    ,8HDEBUG   /
      DATA UHF,ROHF/8HUHF     ,8HROHF    /
C
C     ----- THIS IS THE MAIN 2E- GRADIENT DRIVER -----
C
      DBG = EXETYP.EQ.DEBUG
      OUT = EXETYP.EQ.GRD2.OR.NPRINT.EQ.-4
      SOME = MASWRK  .AND.  NPRINT.NE.-5
      IF(SOME) WRITE (IW,9008)
C
C        INITIALIZE PARALLEL
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT = -1
      MINE = -1
C
C     ----- SET STARTING PARAMETERS -----
C
      HONDO = .TRUE.
      POPLE = igrdtyp.ne.2
C
C      CUTOFF IS THE SCHWARZ SCREENING CUT OFF
C      DABCUT IS THE TWO PARTICLE DENSITY CUT OFF
C
      CUTOFF=TENM9
      IF(.NOT.POPLE) CUTOFF=CUTOFF/TEN
      CUTOFF2=CUTOFF/2.0D+00
C
      ZBIG = 0.0D+00
      DO ISH=1,NSHELL
         I1=KSTART(ISH)
         I2=I1+KNG(ISH)-1
         DO IG=I1,I2
            IF(EX(IG).GT.ZBIG) ZBIG = EX(IG)
         ENDDO
      ENDDO
      DABCUT=TENM11
      IF(ZBIG.GT.1.0D+06) DABCUT = DABCUT/TEN
      IF(ZBIG.GT.1.0D+07) DABCUT = DABCUT/TEN
C
C      VTOLS ARE CUT OFFS USED BY THE POPLE PACKAGE
C      CURRENT VALUES ARE FROM HONDO 8, SEE G92 FOR OTHER POSSIBILITIES
C
      VTOL1 = TENM12
      VTOL2 = TENM12
      VTOLS = TENM20
      DTOL = TEN**(-ITOL)
      RTOL = RLN10*ITOL
C
C      INITIALIZE THE INTEGRAL BLOCK COUNTERS TO ZERO
C
      IISKIP1= 0
      IISKIP2= 0
      IDID = 0
C
      UHFTYP=SCFTYP.EQ.UHF .OR. SCFTYP.EQ.ROHF
C
C     ----- SET POINTERS FOR PARTITIONING OF MEMORY -----
C
c     L1 = NUM
      L2 = (NUM*NUM+NUM)/2
      NSH2=(NSHELL*NSHELL+NSHELL)/2
C
      DO 100 I = 1,NUM
         IA(I) = (I*I-I)/2
  100 CONTINUE
C
C     ----- READ IN 1E-GRADIENT -----
C
c     CALL DAREAD(IDAF,IODA,DE,3*NAT,3,0)
c     IF (GOPARR) CALL DSCAL(3*NAT,ONE/NPROC,DE,1)
c     call vclr(DE,1,3*NAT)
c     DE is set in ESDNRDER
c
C              CALCULATE THE LARGEST SHELL TYPE
C
      CALL BASCHK(MAXTYP)
      MAXSHL = LENSHL(MAXTYP+1)
C              DO AT LEAST AN L SHELL
      IF (MAXSHL.LT.4) MAXSHL=4
C
C       IF WE ARE USING THE POPLE PACKAGE AND DO NOT HAVE ANY SHELLS
C       LARGER THAN AN L-SHELL THEN SKIP THE SETUP FOR THE RYS PACKAGE
C
      IF (POPLE.AND.MAXTYP.LT.2) HONDO = .FALSE.
C
C     FIGURE OUT THE MEMORY WE NEED FOR STORING DENSITY MATRIX
C     AND OTHER WAVEFUNCTION INFORMATION. -JKDMEM- ALLOCATES
C     MEMORY FOR DERIVATIVE COMPUTATION AND 2ND ORDER DENSITY
C     AFTER -LENGTH- WORDS.
C
      LENGTH=L2
      IF (NTMO.GT.0) LENGTH = LENGTH + L2
C
C       CALCULATE THE AMOUNT OF MEMORY NEEDED AND SET THE POINTERS
C       FOR BOTH PACKAGES
C
      CALL VALFM(LOADFM)
      CALL JKDMEM(1,LOADFM,IADDR,LENGTH,MINXYZ,MAXXYZ,MINVEC,POPLE,
     *            .false.)
c     CALL DDI_SYNC(1147)
C
C     ----- CARRY OUT SET UP TASKS -----
C
      NEED=IADDR-LOADFM
      CALL GETFM(NEED)
      IF (EXETYP .EQ. CHECK) GO TO 600
C
C     ----- READ WAVEFUNCTION INFORMATION -----
C
c     IF (NTMO.GT.0) THEN
c        CALL DAREAD(IDAF,IODA,X(LVEC),L2,79,0)
c     END IF
c     CALL DDI_SYNC(1148)
C
C        READ IN THE EXCHANGE INTEGRALS FROM DISK. IF THEY WERE NOT
C        PREVIOUSLY COMPUTED, THEN JUST SET THE ARRAY TO ONE, WHICH
C        EFFECTIVELY DEACTIVATES THE SCHWARZ SCREENING
C
      IF(ISCHWZ.EQ.1) THEN
         CALL DAREAD(IDAF,IODA,X(IXCH),NSH2,54,0)
      ELSE
         DO 400 I=0,NSH2-1
            X(IXCH+I) = ONE
  400    CONTINUE
      END IF
C
C     ----- GET SYMMETRY MAPPING OF SHELLS -----
C
c     JKDSET does absolutely nothing useful for us?
c     CALL JKDSET
C
C        SET UP THE 1-ELECTRON CHARGE DISTRIBUTION
C
      IF (HONDO) CALL OEDHND(X(INIJG),X(ICHRG))
C
C        SQUARE DTOL FOR USE IN JKDSPD
C
      DTOL = DTOL*DTOL
      NC=1
      LDF=1
C
C     ----- I SHELL -----
C
      DO 560 II = ncursh+1,NSHELL
C
C     ----- J SHELL -----
C
        DO 540 JJ = ncursh+1,II
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 540
          END IF
C
C        GET IJ CHARGE DISTRIBUTION
C        REALLY WE JUST SET THE POINTERS TO THE CHARGE DISTRIBUTION
C
          IF (HONDO) THEN
            IIJJ=IA(MAX0(II,JJ))+MIN0(II,JJ)
            CALL OEDRD(X(INIJG),NIJ,NIJ0,IIJJ)
            IF(NIJ.EQ.0) GO TO 540
          END IF
C
C     ----- K SHELL -----
C
          DO 520 KK = 1,ncursh 
C
C     ----- L SHELL -----
C
            DO 500 LL = 1,KK
C
C     ----- GO PARALLEL! -----
C
              IF ((.NOT.NXT) .AND. GOPARR) THEN
                 IPCOUNT = IPCOUNT + 1
                 IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 500
              END IF
C
C     ----- CALCULATE Q4 FACTOR FOR THIS GROUP OF SHELLS -----
C
              Q4 = 1
C
C     ----- DECIDE ON DERIVATIVE INTEGRAL METHOD -----
C     ANY PURELY SP SET OF SHELLS CAN BE DONE WITH THE FASTER
C     POPLE/SCHLEGEL ROTATION ALGORITHM.  INTEGRALS INVOLVING
C     D AND HIGHER FUNCTIONS MUST USE RYS POLYNOMIAL CODE.
C
              POPLE = .TRUE.
              IF(igrdtyp.EQ.2)   POPLE=.FALSE.
              IF(KTYPE(II).GT.2) POPLE=.FALSE.
              IF(KTYPE(JJ).GT.2) POPLE=.FALSE.
              IF(KTYPE(KK).GT.2) POPLE=.FALSE.
              IF(KTYPE(LL).GT.2) POPLE=.FALSE.
C
C         IMPLEMENT INTEGRAL SCREENING HERE USING EXCHANGE INTEGRALS
C         SEE H.HORN, H.WEISS, M.HAESER, M.EHRIG, R.AHLRICHS
C             J.COMPUT.CHEM. 12, 1058-1064(1991) 
C         REGARDING THE ESTIMATION FORMULA (31) THAT IS USED HERE.
C
              IJIJ=IA(MAX0(II,JJ))+MIN0(II,JJ)
              KLKL=IA(MAX0(KK,LL))+MIN0(KK,LL)
              GMAX=(X(IXCH+IJIJ-1)*X(IXCH+KLKL-1))
C
C                COARSE SCREENING, ON JUST THE INTEGRAL VALUE
C
              IF (GMAX.LT.CUTOFF) THEN
                 IISKIP1 = IISKIP1+1
                 GO TO 500
              END IF
C
           ISH=II
           JSH=JJ
           KSH=KK
           LSH=LL
C
           IF (POPLE) THEN
             AX1=PT5
             IF(ISH.NE.JSH) AX1=AX1+AX1
             IF(KSH.NE.LSH) AX1=AX1+AX1
c            IF(ISH.NE.KSH.OR.JSH.NE.LSH) AX1=AX1+AX1
             AX1=AX1+AX1
             Q4 = Q4*AX1
             INEW=ISH
             JNEW=JSH
             KNEW=KSH
             LNEW=LSH
             IMAX=KTYPE(INEW)-1
             JMAX=KTYPE(JNEW)-1
             KKKMAX=KTYPE(KNEW)-1
             LMAX=KTYPE(LNEW)-1
             IF (IMAX.LT.JMAX) THEN
               INEW = JSH
               JNEW = ISH
             END IF
             IF (KKKMAX.LT.LMAX) THEN
               KNEW = LSH
               LNEW = KSH
             END IF
             IF ((IMAX+JMAX).LT.(KKKMAX+LMAX)) THEN
               ID = INEW
               INEW = KNEW
               KNEW = ID
               ID = JNEW
               JNEW = LNEW
               LNEW = ID
             END IF
             IMAX=3*(KTYPE(INEW)-1)+1
             JMAX=3*(KTYPE(JNEW)-1)+1
             KKKMAX=3*(KTYPE(KNEW)-1)+1
             LMAX=3*(KTYPE(LNEW)-1)+1
             JTYPE=(IMAX+JMAX+KKKMAX+KKKMAX+LMAX-2)/3
             IAT = KATOM(INEW)
             JAT = KATOM(JNEW)
             KAT = KATOM(KNEW)
             LAT = KATOM(LNEW)
             IF ((IAT.EQ.JAT).AND.(IAT.EQ.KAT).AND.(IAT.EQ.LAT))
     1           GO TO 500
           ELSE
C
C     ----- GET -KL- CHARGE DISTRIBUTION -----
C       ACTUALLY JUST THE POINTERS
C
              KKLL=IA(MAX0(KK,LL))+MIN0(KK,LL)
              CALL OEDRD(X(INIJG),NKL,NKL0,KKLL)
              IF(NKL.EQ.0) GO TO 500
C
C     ----- SELECT CENTERS FOR DERIVATIVES -----
C
              CALL JKDATM(ISH,JSH,KSH,LSH)
              IF(SKIPI.AND.SKIPJ.AND.SKIPK.AND.SKIPL) GO TO 500
C
C     ----- SET INDICES FOR SHELL BLOCK -----
C
              CALL JKDSHL(ISH,JSH,KSH,LSH)
              CALL JKDNDX(X(IIJKLG))
              INEW = ISH
              JNEW = JSH
              KNEW = KSH
              LNEW = LSH
           END IF
C
C     ----- OBTAIN 2 BODY DENSITY FOR THIS SHELL BLOCK -----
C
c        IF (NTMO.GT.0) THEN
c           CALL DABPAU(INEW,JNEW,KNEW,LNEW,UHFTYP,X(LDEN),
c    *                  X(LVEC),X(IDAB),DABMAX,Q4,POPLE)
c        else
            CALL DABCLU(INEW,JNEW,KNEW,LNEW,UHFTYP,DJ,
     *                  DI,X(IDAB),DABMAX,Q4,POPLE,l1i)
c        END IF
c        DA and DB are flipped around, since the first 2 indices are for J
c        and the other 2 for I. This agrees with FMO2ei. 
C
C                FINE SCREENING, ON INTEGRAL VALUE TIMES DENSITY FACTOR
C
         IF(DABMAX*GMAX.LT.CUTOFF2) THEN
            IISKIP2 = IISKIP2+1
            GO TO 500
         END IF
C
C     ----- EVALUATE DERIVATIVE INTEGRAL, AND ADD TO THE GRADIENT -----
C
         IDID = IDID+1
         IF(POPLE) THEN
            CALL JKDG80(DABMAX,INEW,JNEW,KNEW,LNEW,
     *                  JTYPE,IAT,JAT,KAT,LAT)
         ELSE
            CALL JKDSPD(NIJ0,NKL,NKL0,X(ICHRG),
     *      X(IGINT),X(IFINT),X(ISINT),X(IIJKLG),X(IGIJKL),
     1      X(IGNKL),X(IGNM),X(IXY),X(IXZ),X(IYZ),X(IX),X(IY),X(IZ),
     2      X(ISJ),X(ISK),X(ISL),X(IB00),X(IB01),X(IB10),X(IC00),
     3      X(ID00),X(IF00),X(IDIJ),X(IDKL),X(IDIJSI),X(IDIJSJ),
     4      X(IDKLSK),X(IDKLSL),X(IABV),X(ICV),X(IRW),X(IAAI),X(IAAJ),
     5      X(IBBK),X(IBBL),X(IFI),X(IFJ),X(IFK),X(IFL),X(ISII),X(ISJJ),
     6      X(ISKK),X(ISLL),X(ISIJ),X(ISIK),X(ISIL),X(ISJK),X(ISJL),
     7      X(ISKL),X(IDAB),MAXXYZ,FC,NC,DF,LDF,NBF,DDA,Q4,MINVEC,
     8      DABCUT,DABMAX)
         END IF
C
C     ----- END OF *SHELL* LOOPS -----
C
  500 CONTINUE
  520 CONTINUE
  540 CONTINUE
      IF (TIM .GE. TIMLIM) GO TO 600
  560 CONTINUE
C
C     ----- FINISH UP THE FINAL GRADIENT -----
C
C           GLOBAL SUM CONTRIBUTIONS FROM EACH NODE
C
      IF (GOPARR) THEN
         IF (NXT) CALL DDI_DLBRESET
c        CALL DDI_GSUMF(1600,DE,3*NAT)
c        Do not global sum: it will be done elsewhere (FMOX).
c        CALL DDI_GSUMI(1601,IISKIP1,1)
c        CALL DDI_GSUMI(1601,IISKIP2,1)
c        CALL DDI_GSUMI(1602,IDID   ,1)
         idum=0
         call DDI_nsumi(1601,IISKIP1,IISKIP2,IDID,idum,3)
      END IF
C
C           SYMMETRIZE SKELETON GRADIENT VECTOR
C
c     CALL SYMEG(DE)
c     Do not symmetrise: no symmetry. If one is to symmetrise, other 
c     contributions fron TVDER, SDER and ESDNRDER are to be pre-added.
C
C     ----- DEALLOCATE MEMORY -----
C
  600 CONTINUE
c     IF(LAST.GT.0) LAST=0
      CALL RETFM(NEED)
      if(nprint.ne.-5) then
        IF (MASWRK) THEN
         IF(SOME) WRITE(IW,9999) IISKIP1,IISKIP2,IDID
         WRITE(IW,FMT='('' ...... END OF 2-ELECTRON GRADIENT ......'')')
        END IF
        CALL TEXIT(1,4)
      endif
c     WRITE(IW,9088)
c     CALL EGOUT(DE,NAT)
      RETURN
C
 9008 FORMAT(/10X,22(1H-)/10X,22HGRADIENT OF THE ENERGY/10X,22(1H-))
c9088 FORMAT(/10X,36(1H-)/10X,'GRADIENTincluding 2e terms',/10X,36(1H-))
 9999 FORMAT(1X,'THE COARSE/FINE SCHWARZ SCREENINGS SKIPPED ',I10,'/',
     *          I10,' BLOCKS.'/
     *     1X,'THE NUMBER OF GRADIENT INTEGRAL BLOCKS COMPUTED WAS',I10)
c--      END
c--C*MODULE fmoint  *DECK mapdist
c--      SUBROUTINE mapdist(map,indat,indatg,iaglob,ifg,jfg,kfg,lfg)
c--      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
c--      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
c--      logical IANDJ,bimer(3)
c--      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
c--     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
c--      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
c--     *                CF(MXGTOT),CG(MXGTOT),
c--     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
c--     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
c--      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
c--      dimension map(*),indat(*),indatg(natfmo,*),iaglob(*)
c--      data bimer/3*.false./
c--c
c--      if(jfg.ne.0) then
c--        bimer(1)=fmodist(ifg,0,0,jfg).eq.0
c--        if(kfg.ne.0) then
c--          bimer(2)=fmodist(ifg,0,0,kfg).eq.0
c--          bimer(3)=fmodist(jfg,0,0,kfg).eq.0
c--          if(bimer(1).and.bimer(2)) bimer(3)=.true.
c--          if(bimer(1).and.bimer(3)) bimer(2)=.true.
c--          if(bimer(2).and.bimer(3)) bimer(1)=.true.
c--        endif
c--      endif
c--      DO 720 II = 1,NSHELL
c--         I = KATOM(II)
c--         MINI = KMIN(II)
c--         MAXI = KMAX(II)
c--         LOCI = KLOC(II)-MINI
c--         iatg=iaglob(I)
c--         DO 700 JJ = 1,II
c--            J = KATOM(JJ)
c--            MINJ = KMIN(JJ)
c--            MAXJ = KMAX(JJ)
c--            LOCJ = KLOC(JJ)-MINJ
c--            jatg=iaglob(j)
c--            IANDJ = II .EQ. JJ
c--            MAX = MAXJ
c--            DO 620 I = MINI,MAXI
c--               LI = LOCI+I
c--               IN = (LI*(LI-1))/2
c--               IF (IANDJ) MAX = I
c--               DO 600 J = MINJ,MAX
c--                  LJ = LOCJ+J
c--                  JN = LJ+IN
c--                  map(jn)=int(1.0D+02*fmossdist(indat(iatg),indat(jatg),
c--     *             indatg(iatg,1),indatg(jatg,1),ifg,jfg,kfg,lfg,bimer))
c--  600          CONTINUE
c--  620       CONTINUE
c--  700    CONTINUE
c--  720 CONTINUE
c--      loop=0
c--      do i=1,num
c--        write(6,9000) i,(map(loop+j),j=1,i)  
c--        loop=loop+i 
c--      enddo
c-- 9000 format(9999I4,/)
c--      RETURN
         END
C*MODULE fmoint  *DECK indsalc
      SUBROUTINE indsalc(salc,ind,nsalc,l1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXSH=5000,MXGTOT=20000,MXATM=2000)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      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
      dimension salc(num,num),ind(num)
c
c     Index the Z (SALC) matrix. ind(i) contains the number of the first 
c     non-zero element in the column i (which is the AO index). 
c     SALC matrix is sorted as s,p,d etc, so that columns are interchanged
c     except for some trivial cases.
c
c     if(iand(ixesp,128).eq.0.or.nbdfg.eq.0) return
      call viclr(ind,1,l1)
c     loop over vectors (columns)
      do lj=1,nsalc
c       loop over rows
        DO II=1,NSHELL
          MINI=KMIN(II)
          MAXI=KMAX(II)
          LOCI=KLOC(II)-MINI
          do i=mini,maxi
            LI = LOCI+i
            if(salc(li,lj).ne.0) then
c             write(6,*) 'Found ',li,lj
              ind(lj)=li
              goto 100
            endif
          enddo
        enddo
  100   continue
      enddo
      return
      end
C*MODULE fmoint  *DECK copyZbl
      SUBROUTINE copyZbl(nao,nmo,l1,iloc,ind,Zbl,ndimZbl,Z,ndimZ,nsalc)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension ind(*),Zbl(ndimZbl,nmo),Z(ndimZ,nao)
c
c     Copy an atomic block into the Z (SALC) matrix. The block is obtained
c     from atomic calculations of diagonalising projection operators and the 
c     purpose of the block is to remove some basis functions from the
c     Z matrix, very much like removing spherical contaminants. 
c 
      icopied=0
      do i=1,l1
c       SALC matrix is sorted according to L values, so find the proper place 
        if(ind(i).ge.iloc.and.ind(i).lt.iloc+nao) then
          icopied=icopied+1
          call dcopy(nao,Zbl(1,icopied),1,Z(iloc,i),1)
          if(icopied.gt.nmo) then
c           This means erasing an AO from Z (=SALC).
c           ind(i)=0 means zero vector in Z, to be shifted back later
            ind(i)=0
            nsalc=nsalc-1
          endif
        endif
      enddo
      return
      end
C*MODULE fmoint  *DECK pushZback
      SUBROUTINE pushZback(Z,ind,indnew,l1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension Z(l1,l1),ind(l1),indnew(l1)
c
c     Sort the Z (SALC) matrix by pushing all zero vectors toward the end
c
      iloop=0
c     first renumber nonzero indices
      do i=1,l1
        if(ind(i).ne.0) then
          iloop=iloop+1
          indnew(iloop)=i
        endif
      enddo 
c     next, renumber zero indices
      do i=1,l1
        if(ind(i).eq.0) then
          iloop=iloop+1
          indnew(iloop)=i
        endif
      enddo
      call REORDR(Z,indnew,l1,l1)
c     write(6,*) 'final Z'
c     call prsq(Z,l1,l1,l1)
      return
      end
C*MODULE fmoint  *DECK FMOqmt
      SUBROUTINE FMOqmt(s,q,v,e,p,scr,wrk,L1,L2,L3,nzapped,nqmt) 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      dimension s(l2),q(l1,l1),v(l1,l1),e(l1),p(l2),scr(l1,8),wrk(l1)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      common /fmorun/ espscf,IDAFMO,icurlay,icurfg,jcurfg,kcurfg,ncursh,
     *                ngau,icurpop,ifmostp,moncor,nat1e,needr,modrst,
     *                norbproj
c
c     Construct the new Z (=SALC) matrix for FMO, excluding
c     a) spherical contaminants (same as for regular (not FMO) runs.
c     b) some basis functions corresponding to MOs along the fractioned bonds.
c     If none functions are zapped according to b) then just generated
c     regular Q matrix and quit.
c
c     On the input: P, the projection operator matrix, S overlap integrals.
c     On the output: V, the new SALC matrix of size (l1,nqmt).
c
c     1. find orthogonal basis Qt * S * Q =1
c     To force QMTSYM to generate the original Q matrix, reset ifmostp.
c
      ifmosav=ifmostp
      ifmostp=1
      if(nzapped.eq.0) then
c       put the results into V.
        call QMTSYM(S,Q,v,E,SCR,WRK,L0,L1,L2,L3,DBG)
        nqmt=l0
        ifmostp=ifmosav
        return
      else
        call QMTSYM(S,v,Q,E,SCR,WRK,L0,L1,L2,L3,DBG)
      endif
      ifmostp=ifmosav
c
c     2. transform the projection operator P to the new basis (=P')
c
      call TFTRI(S,p,Q,WRK,l0,l1,l1)
c
c     3. diagonalise P'. nzapped is the number of eigenvectors to be removed.
c     They always come last because their eigenvalues are nearly equal to
c     B, which is 1e+6 or so and other eigenvalues are nearly 0, since
c     P and P' have an approximate rank of nzapped. 
c
      CALL GLDIAG(L1,L0,l0,s,SCR,E,v,IERR,WRK)
      IF (IERR .NE. 0) CALL ABRT
c
c     4. transform the chosen eigenvectors back to AO basis 
c     and use the eigenvectors as the new SALC matrix!
c     (clear the removed part, just in case).
c
      CALL TFSQB(v,Q,wrk,l0,l1,l1)
      nqmt=l0-nzapped
      write(6,*) 'wwwQ',nqmt,l0,l1
      call prsq(v,l1,l1,l1)
      call vclr(v(1,nqmt+1),1,nzapped*l1)
c
c     This produced the new SALC matrix which happens to be also a Q matrix,
c     since it satisfies Qt * S * Q = 1. 
c
      if(maswrk) WRITE(IW,9003) nzapped 
      CALL DAREAD(IDAF,IODA,s,l2,12,0)
      call TFTRI(p,s,V,WRK,nqmt,l1,l1)
      call prtri(p,nqmt)
c
      return
 9003 FORMAT(1X,'NUMBER OF FMO redundant MOS DROPPED=     ',I6)
      end
C*MODULE fmoint  *DECK convorb
      subroutine convorb(ifg,l1,da,l1c,orbconv,ibfconv,mapconv,iaglob,
     *                   izbas)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000,MXAO=8192,MAXL=5,MaxLay=5,MXSFMO=MaxLay*40,
     *           MXGFMO=MaxLay*100,MXAFMO=MaxLay*10)
      dimension da(*),orbconv(*),ibfconv(2,*),mapconv(*),
     *          locat(MXATM),locmap(MXAFMO),iaglob(*),izbas(*)
      logical GOPARR,DSKWRK,MASWRK
      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,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      common /fmoshl/ natl,numl,fzan(MXAFMO),fC(MXGFMO,MAXL),fEX(MXGFMO)
     *               ,lshell,LSTART(MXSFMO),LATOM(MXSFMO),LTYPE(MXSFMO),
     *                LNG(MXSFMO),LMIN(MXSFMO),LMAX(MXSFMO),llay(MXAFMO)
c
c     Convert density for one fragment from one basis set to another.
c     Arguments:
c     ibfconv(1,nat in $data) the number of old basis function fo reach atom
c     ibfconv(2,nat in $data) the number of new basis function fo reach atom
c     mapconv(1) array givinig mapping of the new basis in terms of the old,
c     e.g. H2O, 631G(old) -> 631G*(new) (provided in $DATA the order is O, H):
c     ibfconv(1)=9,15, 1,1
c     mapconv(1)=1,2,3,4,5,6,7,8,9,0,0,0,0,0,0, 1
c     locat, locmap: temporary arrays.  
c
c     First, find the new l1c
      l1o=0
      l1c=0
      do iat=1,nat
        iz=ian(iat)
        ibas=izbas(iaglob(iat))
        iatd=0
        do jat=1,natl
          jz=int(fzan(jat)+0.5D+00)
          jbas=llay(jat)/100
          if(iz.eq.jz.and.ibas.eq.jbas) then
            iatd=jat
            goto 100
          endif
        enddo
  100   continue
        if(iatd.eq.0) call abrt
        locat(iat)=iatd
        l1o=l1o+ibfconv(1,iatd)
        l1c=l1c+ibfconv(2,iatd)
      enddo
c     Fill locmap
      map0=1
      do iat=1,natl
        locmap(iat)=map0
        map0=locmap(iat)+ibfconv(2,iat)
      enddo
      if(maswrk) write(iw,9000) ifg,l1,l1c 
      if(l1.ne.l1o) then
        if(maswrk) write(iw,*) 'Wrong maps in convorb',l1,l1o
        call abrt
      endif
      l2c=(l1c*l1c+l1c)/2
      call vclr(orbconv,1,l2c)
c
      iloco=0
      ilocn=0
      do iat=1,nat
        iatd=locat(iat)
        m1i=ibfconv(2,iatd)
        imap=locmap(iatd)
        do i=1,m1i
          ilocn=ilocn+1
          ilocm=mapconv(imap+i-1)
c         If i-row in the new basis set maps to something in the old
          if(ilocm.ne.0) then
            iold=iloco+ilocm
            jloco=0
            jlocn=0
            do jat=1,nat
              jatd=locat(jat)
              m1j=ibfconv(2,jatd)
              jmap=locmap(jatd)
              do j=1,m1j
                jlocn=jlocn+1
                jlocm=mapconv(jmap+j-1)
c         If j-column in the new basis set maps to something in the old
                if(jlocm.ne.0) then
                  jold=jloco+jlocm
                  if(ilocn.ge.jlocn) then
c                   It is not safe to use IA for loopn: IA is filled for the
c                   old basis which can be (and usually is) smaller than new.
                    loopn=(ilocn*ilocn-ilocn)/2+jlocn
                    loopo=ia(max(iold,jold))+min(iold,jold)
                    orbconv(loopn)=da(loopo)
                  endif
                endif
              enddo
              jloco=jloco+ibfconv(1,jatd)
            enddo
          endif
        enddo
        iloco=iloco+ibfconv(1,iatd)
      enddo
      write(6,*) 'Old orbitals'
      call prtril(da,l1)
      write(6,*) 'New orbitals'
      call prtri(orbconv,l1c)
      return
 9000 format(/1x,'Converting ifg=',I4,' L1=',I4,' to L1=',I4)
      END
C*MODULE fmo     *DECK bconvrec
      subroutine bconvrec(fmozan,izbas,iabdfg,NDAR30,iodcfmo,l1fmo,
     *                    l1fmoc,IDAcFMO,ibfconv,compl1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAXL=5,MaxLay=5,MXSFMO=MaxLay*40,
     *           MXGFMO=MaxLay*100,MXAFMO=MaxLay*10)
      logical compl1
      COMMON /RAIOLN/ JRECLN(10),JRECST(10)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,convfg,cnvdmp,coroff,IXESP,mxitfg,nguess,
     *                NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c
      common /fmoshl/ natl,numl,fzan(MXAFMO),fC(MXGFMO,MAXL),fEX(MXGFMO)
     *               ,lshell,LSTART(MXSFMO),LATOM(MXSFMO),LTYPE(MXSFMO),
     *                LNG(MXSFMO),LMIN(MXSFMO),LMAX(MXSFMO),llay(MXAFMO)
      dimension fmozan(*),izbas(*),iabdfg(*),iodcfmo(*),ibfconv(2,*)
c
c     Fill in initial (index) record for the orbital conversion.
c
c     if(maxl1.ne.maxl1c) then
      if(compl1) then
        l1bd=0
        do ibdg=1,nbdfg
          ia=abs(iabdfg(ibdg))
          iz=int(fmozan(ia)+0.5D+00)
          ibas=izbas(ia) 
          iatd=0
          do jat=1,natl
            jz=int(fzan(jat)+0.5D+00)
            jbas=llay(jat)/100
            if(iz.eq.jz.and.ibas.eq.jbas) then
              iatd=jat
              goto 100
            endif
          enddo
  100     continue
          if(iatd.eq.0) call abrt
          nao=ibfconv(2,iatd)
          l1bd=l1bd+nao
        enddo
        l1fmoc=l1fmoc-l1bd
      else
        l1bd=0
        l1fmoc=l1fmo
      endif
c
      iodcfmo(NDAR30+1)=modorb
      iodcfmo(NDAR30+2)=l1fmoc
      iodcfmo(NDAR30+3)=irstlay
      iodcfmo(NDAR30+4)=1
c     1 means next iteration after restart is number 1 (converted orbitals
c     are somewhat alike initial guess, iteration 0).
c     Only one set of orbitals is written: possible problem with MP2 restarts.
      iodcfmo(NDAR30+5)=JRECST(IDAcFMO/10)
c     call icopy(nlayer,idmrec0,1,iodcfmo(NDAR30+6),1)
      do i=1,nlayer
        iodcfmo(NDAR30+5+i)=1
      enddo
c     write(6,*) 'Saving for restart:',l1fmoc,l1bd,iodcfmo(NDAR30+5)
      return
      END
c
C*MODULE fmoint  *DECK fmoprc
      SUBROUTINE fmoprc(ida,iblock,l2,da,dgrid,minx,NXG,miny,NYG,minz,
     *                  NZG,ixmin,ixmax,iymin,iymax,izmin,izmax,jxmin,
     *                  jxmax,jymin,jymax,jzmin,jzmax,ORIGIN,UX,UY,UZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION da(l2),dgrid(minz:NZG,miny:NYG,minx:NXG),
     *          ORIGIN(3),UX(3),UY(3),UZ(3)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,IANDJ,NORM,DOUBLE,nxt,testyz,testz2
C
      DIMENSION IX(35),IY(35),IZ(35),JX(35),JY(35),JZ(35)
      DIMENSION DIJ(225),IJX(225),IJY(225),IJZ(225)
      DIMENSION XIN(25),YIN(25),ZIN(25)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),ch(mxgtot),ci(mxgtot),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PRPINT/ XINT0,XINT1,XINT2,XINT3,
     *                YINT0,YINT1,YINT2,YINT3,
     *                ZINT0,ZINT1,ZINT2,ZINT3
      COMMON /XYZORB/ T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ,NM
      COMMON /XYZPRP/ XP,YP,ZP
     *               ,DMX,DMY,DMZ
     *               ,QXX,QYY,QZZ,QXY,QXZ,QYZ
     *               ,QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ
     *               ,OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ
     *               ,OXZZ,OYZZ,OZZZ,OXYZ
     *               ,OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY
     *               ,OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
C
      PARAMETER (SQRT3=1.73205080756888D+00)
      PARAMETER (SQRT5=2.23606797749979D+00)
      PARAMETER (SQRT7=2.64575131106459D+00)
      PARAMETER (ZERO=0.0D+00, tiny=1.0D-08)
      PARAMETER (ONE=1.0D+00, two=2.0D+00)
      PARAMETER (RLN10=2.30258D+00)
C
      DATA JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     *          3, 0, 0, 2, 2, 1, 0, 1, 0, 1,
     *          4, 0, 0, 3, 3, 1, 0, 1, 0, 2,
     *          2, 0, 2, 1, 1/
      DATA IX / 1, 6, 1, 1,11, 1, 1, 6, 6, 1,
     *         16, 1, 1,11,11, 6, 1, 6, 1, 6,
     *         21, 1, 1,16,16, 6, 1, 6, 1,11,
     *         11, 1,11, 6, 6/
      DATA JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     *          0, 3, 0, 1, 0, 2, 2, 0, 1, 1,
     *          0, 4, 0, 1, 0, 3, 3, 0, 1, 2,
     *          0, 2, 1, 2, 1/
      DATA IY / 1, 1, 6, 1, 1,11, 1, 6, 1, 6,
     *          1,16, 1, 6, 1,11,11, 1, 6, 6,
     *          1,21, 1, 6, 1,16,16, 1, 6,11,
     *          1,11, 6,11, 6/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     *          0, 0, 3, 0, 1, 0, 1, 2, 2, 1,
     *          0, 0, 4, 0, 1, 0, 1, 3, 3, 0,
     *          2, 2, 1, 1, 2/
      DATA IZ / 1, 1, 1, 6, 1, 1,11, 1, 6, 6,
     *          1, 1,16, 1, 6, 1, 6,11,11, 6,
     *          1, 1,21, 1, 6, 1, 6,16,16, 1,
     *         11,11, 6, 6,11/
C
C       CALCULATE INTEGRALS NECESSARY FOR EVALUATION OF PROPERTIES.
c       This is a clone of PRCALC.
c
c     Accumulate the electron density for the total system on the grid
c     determined by ixmin,ixmax,iymin,iymax,izmin,izmax and excluding (iblock>1)
c     the part given by jxmin,jxmax,jymin,jymax,jzmin,jzmax.
c
c     Warning: this subroutine uses a fixed value of ITOL. It is thought 
c     (and checked) that the total density on the grid if is used mainly
c     for plotting does not need extra fine accuracy while the expenses
c     are considerate. 
c
c     Written by D. G. Fedorov using suggestions by Y. Inadomi.
c     parstat: GroupFull.
c     
      if(ida.eq.0.or.ixmax.lt.ixmin.or.iymax.lt.iymin.or.izmax.lt.izmin)
     *  return
c
      if(maswrk) write(iw,9000)ixmax-ixmin+1,iymax-iymin+1,izmax-izmin+1
C
c     itol0=20
      itol0=18
c     cutoff=1.0D-06
c     cutoff is the AO density threshold (see below for more explanation).
c     The total electron density at a grid point is the sum of { the AO density 
c     times the product of the values of two AOs at the point }. 
c     Since the latter product is (usually much) smaller than 1, cutoff clearly
c     governs the overall accuracy. Cube files print only 5 digits so at most
c     one wants 5 digits of accuracy (in the relative sense). Given the sum, it
c     should be sufficient to use a cutoff of slightly smaller than that. 
c     Note that for exponents (itol0) fairly conservative value must be used as
c     its effect is less controllable.
c     Preliminary tests sgowed that cutoff=1e-6 has reasonable accuracy
c     (at worst one in the last printed digit) but almost no savings.
c
      TOL = RLN10*itol0
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      NXT = IBTYP.EQ.1
      NEXT  = -1
      kount = -1
      testz2=abs(UZ(1)).lt.tiny.and.abs(UZ(2)).lt.tiny
      testyz=abs(UY(1)).lt.tiny.and.abs(UZ(1)).lt.tiny
      fda=ida
c
c     Prepare density for easy tracing: double off-diagonal elements.
c
      l1=num
      loop=0
      do i=1,l1
        do j=1,i-1
          loop=loop+1
          da(loop)=da(loop)*two 
        enddo
c       skip the diagonal
        loop=loop+1
      enddo
C
C       LOOP OVER SHELLS II
C
      DO 500  II=1,NSHELL
        I    = KATOM(II)
        XI   = C(1,I)
        YI   = C(2,I)
        ZI   = C(3,I)
        I1   = KSTART(II)
        I2   = I1 + KNG(II) - 1
        LIT  = KTYPE(II)
        MINI = KMIN(II)
        MAXI = KMAX(II)
        LOCI = KLOC(II) - MINI
C
C         LOOP OVER SHELLS JJ
C
        DO 500  JJ=1,II
c     
          if(goparr) then
            kount=kount+1
            if(nxt) then
              IF(kount.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
              if(NEXT.ne.kount) goto 500
            else 
              if(MOD(kount,NPROC).NE.me) goto 500
            endif
          endif
C     
          J    = KATOM(JJ)
          XJ   = C(1,J)
          YJ   = C(2,J)
          ZJ   = C(3,J)
          J1   = KSTART(JJ)
          J2   = J1 + KNG(JJ) - 1
          LJT  = KTYPE(JJ)
          MINJ = KMIN(JJ)
          MAXJ = KMAX(JJ)
          LOCJ = KLOC(JJ) - MINJ
C       
          RR     = (XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2
          IANDJ  = II.EQ.JJ
C       
C           PREPARE INDICES FOR PAIRS OF (I,J) ORBITALS
C       
          IJ = 0
          MAX = MAXJ
          DO 100  I=MINI,MAXI
            NX = IX(I)
            NY = IY(I)
            NZ = IZ(I)
C         
            IF (IANDJ) MAX = I
            DO 100  J=MINJ,MAX
              IJ = IJ+1
              IJX(IJ) = NX+JX(J)
              IJY(IJ) = NY+JY(J)
              IJZ(IJ) = NZ+JZ(J)
  100     CONTINUE
C
c         DO 120  I=1,225
c 120       WINT(I) = ZERO
C
C           LOOP OVER PRIMITIVES IG
C
          JGMAX = J2
          DO 400  IG=I1,I2
            AI  = EX(IG)
            CSI = CS(IG)
            CPI = CP(IG)
            CDI = CD(IG)
            CFI = CF(IG)
            CGI = CG(IG)
C      
C             LOOP OVER PRIMITIVES JG
C    
            IF (IANDJ) JGMAX = IG
            DO 400 JG=J1,JGMAX
              AJ  = EX(JG)
              CSJ = CS(JG)
              CPJ = CP(JG)
              CDJ = CD(JG)
              CFJ = CF(JG)
              CGJ = CG(JG)
C        
              AA  = AI + AJ
              AA1 = ONE/AA
              AX  = (AI*XI + AJ*XJ)*AA1
              AY  = (AI*YI + AJ*YJ)*AA1
              AZ  = (AI*ZI + AJ*ZJ)*AA1
C        
              DUM = AI*AJ*RR*AA1
              IF(DUM .GT. TOL) GO TO 400
              dum0=dum
C
C               CALCULATE DENSITY FACTORS
c               Note that the grid-dependent exponent multiplier FAC
c               was removed from DIJ and used directly much below. 
c               (WINT is proportional to DIJ). 
C    
              DOUBLE = IANDJ.AND.IG.NE.JG
              MAX = MAXJ
              NN  = 0
C             
              DUM1 = ZERO
              DUM2 = ZERO
c             dmax = zero
c             dmax is the max of
c              (AO contraction coefficient product) * (max electron density) 
c             Max is computed for the pair of primitives.
c             In other words, dmax is the max value of the electron density in
c             the primitive basis for all Ml-components with the primitive pair
c             fixed (Ml runs through cartesian components of s,p,d,f,g).
c
              DO 380 I = MINI,MAXI
c               LI = LOCI + I
c               IN = LI*(LI-1)/2
                IF (I.EQ.1) DUM1=CSI
                IF (I.EQ.2) DUM1=CPI
                IF (I.EQ.5) DUM1=CDI
                IF ((I.EQ.8).AND.NORM) DUM1=DUM1*SQRT3
                IF (I.EQ.11) DUM1 = CFI
                IF ((I.EQ.14).AND.NORM) DUM1 = DUM1*SQRT5
                IF ((I.EQ.20).AND.NORM) DUM1 = DUM1*SQRT3
                IF (I.EQ.21) DUM1 = CGI
                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
                IF(IANDJ) MAX = I
                DO 380 J = MINJ,MAX
c                 LJ = LOCJ + J
c                 JN = LJ + IN
                  NN = NN+1
                  IF(J.EQ.1) THEN
                    DUM2 = DUM1*CSJ
                    IF(DOUBLE .AND. I.EQ.1) DUM2 = DUM2 + DUM2
                    IF(DOUBLE .AND. I.GT.1) DUM2 =DUM2+CSI*CPJ
C             
                  ELSE IF(J.EQ.2) THEN
                    DUM2 = DUM1*CPJ
                    IF(DOUBLE) DUM2 = DUM2 + DUM2
C             
                  ELSE IF(J.EQ.5) THEN
                    DUM2 = DUM1*CDJ
                    IF(DOUBLE) DUM2 = DUM2 + DUM2
C             
                  ELSE IF((J.EQ.8).AND.NORM) THEN
                    DUM2 = DUM2*SQRT3
C             
                  ELSE IF (J.EQ.11) THEN
                    DUM2 = DUM1*CFJ
                    IF (DOUBLE) DUM2 = DUM2+DUM2
C             
                  ELSE IF ((J.EQ.14).AND.NORM) THEN
                    DUM2 = DUM2*SQRT5
C             
                  ELSE IF ((J.EQ.20).AND.NORM) THEN
                    DUM2 = DUM2*SQRT3
C             
                  ELSE IF (J.EQ.21) THEN
                    DUM2 = DUM1*CGJ
                    IF (DOUBLE) DUM2 = DUM2+DUM2
C             
                  ELSE IF ((J.EQ.24).AND.NORM) THEN
                    DUM2 = DUM2*SQRT7
C             
                  ELSE IF ((J.EQ.30).AND.NORM) THEN
                    DUM2 = DUM2*SQRT5/SQRT3
C             
                  ELSE IF ((J.EQ.33).AND.NORM) THEN
                    DUM2 = DUM2*SQRT3
C             
                  END IF
                  DIJ(NN) = DUM2
c                 dpij=abs(dum2*da(jn))
c                 if(dpij.gt.dmax) dmax=dpij
  380         CONTINUE
c
c             It is thought better not to include FDA in the cutoff test.
c             The reason is that physically all contributions come as 1.
c             FDA other than 1 come from having monomers subtracted
c             from dimers, so for each FDA>1 with a plus sign there will
c             be about equal contribution subtracted (for 2-body terms
c             and higher, 1-body terms simply add up).  
c             It seems better not to multiply dmax by exp(-dum0) as it will
c             lead to uncontrolled errors due to large multiplicative factors
c             coming from (xi-xp)*(yi-yp)*(zi-zp)*(xj-xp)*(yj-yp)*(zj-zp) 
c             products (computed in DELINT) and ignored in such a test. 
c             One might find the upper bound pmax for them, then 
c             dmax*exp(-dum0)*pmax would have been a great test. 
c             Finding pmax is not as hard as it seems as one knows the
c             largest(smallest) power for each multiplier in the product so
c             there is no need to actually compute all possible products and 
c             find the maximum. The problem is that pmax would depend upon 
c             mx,my and mz and thus will have to be computed inside the grid 
c             loops. Still thinking about 9-fold nested loops here it may be
c             worth doing. 
c
c             IF(dmax.lt.cutoff) GO TO 400
c
              do 490 mx=ixmin,ixmax
                mx0=mx-1
c
c               Extra screening is only possible when y,z-orths are "good".
c               This is perhaps a common situation (i.e., ux,uy and uz are
c               x,y and z-axes, respectively). But it is NOT ALWAYS SO, 
c               so the "goodness" is used only when appropriate. 
c
c               if(abs(UY(1)).lt.tiny.and.abs(UZ(1)).lt.tiny) then
                if(testyz) then
                  XP = ORIGIN(1) + mX0*UX(1)
                  DUM=dum0 + AA*((XP-AX)**2)
                  IF(DUM .GT. TOL) GO TO 490
                endif
                do 480 my=iymin,iymax
                  my0=my-1 
c
c                 Extra screening is only possible when z-orth is "good".
c
c                 if(abs(UZ(1)).lt.tiny.and.abs(UZ(2)).lt.tiny) then
                  if(testz2) then
                    XP = ORIGIN(1) + mX0*UX(1) + mY0*UY(1)
                    YP = ORIGIN(2) + mX0*UX(2) + mY0*UY(2)
                    DUM=dum0 + AA*((XP-AX)**2 + (YP-AY)**2)
                    IF(DUM .GT. TOL) GO TO 480
                  endif
                  do 470 mz=izmin,izmax
c        
c                   Exclude already computed redundant points in a dimer.
c    
                    if(iblock.gt.1 .and. 
     *                     mx.ge.jxmin.and.mx.le.jxmax. and.
     *                     my.ge.jymin.and.my.le.jymax. and.
     *                     mz.ge.jzmin.and.mz.le.jzmax) goto 470 
c
c                   General screening for any orths.
c
                    mz0=mz-1
                    XP = ORIGIN(1) + mX0*UX(1) + mY0*UY(1) + mZ0*UZ(1)
                    YP = ORIGIN(2) + mX0*UX(2) + mY0*UY(2) + mZ0*UZ(2)
                    ZP = ORIGIN(3) + mX0*UX(3) + mY0*UY(3) + mZ0*UZ(3)
                    DUM=dum0 + AA*((XP-AX)**2 + (YP-AY)**2 + (ZP-AZ)**2)
                    IF(DUM .GT. TOL) GO TO 470
C        
                    FAC = EXP(-DUM)
c
c                   fac2= fac*two
C        
c                   CALL INTDEN(LIT,LJT,IJ,IJX,IJY,IJZ,DIJ,WINT)
c
                    IN = -5
                    DO I=1,LIT
                      IN = IN+5
                      NI = I
                      DO J=1,LJT
                        JN = IN+J
                        NJ = J
                        CALL DELINT
                        XIN(JN) = XINT1
                        YIN(JN) = YINT1
                        ZIN(JN) = ZINT1
                      enddo
                    enddo
C   
C                   Contract density with integrals. 
c                   (off-diagonal elements count twice, built in density).
c   
c                   EDENS = TRACEP(da,delta,L1)
c   
c                   dgrid(iz,iy,ix)=dgrid(iz,iy,ix)+ida*EDENS 
C   
                    MAX = MAXJ
                    NN  = 0
                    sum = zero
c                   Tr(Delta*D) (contribution from the pair of primitives).
                    DO 450  I=MINI,MAXI
                      LI = LOCI + I
                      IN = LI*(LI-1)/2
                      IF (IANDJ) MAX = I
                      DO 440  J=MINJ,MAX
                        LJ = LOCJ + J
                        JN = LJ + IN
                        NN = NN+1
                        NX = IJX(NN)
                        NY = IJY(NN)
                        NZ = IJZ(NN)
                        sum=sum+da(jn)*DIJ(NN)*XIN(NX)*YIN(NY)*ZIN(NZ)
c
c                       if(li.eq.lj) then
c                         facd=fac
c                       else
c                         facd=fac2
c                       endif
c                       s=s+facd*da(jn)*DIJ(NN)*XIN(NX)*YIN(NY)*ZIN(NZ)
c                       WINT(I) = WINT(I) + DIJ(I)*XIN(NX)*YIN(NY)*ZIN(NZ)
c                       XVAL(JN) = WINT(NN)
  440                 CONTINUE
  450               CONTINUE
c
                    dgrid(mz,my,mx)=dgrid(mz,my,mx)+fda*fac*sum
c
  470             CONTINUE
  480           CONTINUE
  490         CONTINUE
  400     CONTINUE
C
C        END OF LOOPS OVER PRIMITIVES
C
C
  500 CONTINUE
c
      IF(goparr.and.nxt) CALL DDI_DLBRESET
c
c     Restore density in courtesy
c
      loop=0
      do i=1,l1
        do j=1,i-1
          loop=loop+1 
          da(loop)=da(loop)/two 
        enddo
c       skip the diagonal
        loop=loop+1
      enddo
c
      call timit(1)
c
      RETURN
 9000 format(/1x,'Computing density on the grid ',3I4,/)
      END
C*MODULE fmoint  *DECK nsphel
      function nsphel(ia,ilay)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAXL=5,MAXNZ=137,MaxLay=5,MXSFMO=MaxLay*40,
     *           MXGFMO=MaxLay*100,MXAFMO=MaxLay*10)
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /FMCOM / X(1)
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      common /fmopnt/ lichfg,lmulfg,lidmrec,lfrgnam,llayfrg,lindat,
     *                lnCBS,lfmozan,lfmoc,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(2),lfmode,lnumfrg,lloctat,
     *                liaoglob,lloadm,lfmoge,ldgrid,liodcfmo
      common /fmoshl/ natl,numl,fzan(MXAFMO),fC(MXGFMO,MAXL),fEX(MXGFMO)
     *               ,lshell,LSTART(MXSFMO),LATOM(MXSFMO),LTYPE(MXSFMO),
     *                LNG(MXSFMO),LMIN(MXSFMO),LMAX(MXSFMO),llay(MXAFMO)
c
c     Estimate the number of eliminated spherical contaminants on an atom
c     stored in the basis set library.
c
c     dimension libish(MAXNZ,maxbas,*),libnsh(MAXNZ,maxbas,*)
c     i0=libish(iz,ibas,ilay)
c     n=libnsh(iz,ibas,ilay)
      nsphel=0
      if(ispher.ne.1) return
      ibas=ixftch(x(lizbas),ia)
      iz=int(x(lfmozan+ia-1)+1.0D-03)
      ind=iz+MAXNZ*(ibas-1+maxbas*(ilay-1))
      i0=ixftch(x(llibish),ind)
      n=ixftch(x(llibnsh),ind)
      do i=i0,i0+n-1
c       d starts at 5 (-1), f starts at 11 (-3), g starts at 21 (-6).
        if(lmin(i).eq.5)  nsphel=nsphel+1
        if(lmin(i).eq.11) nsphel=nsphel+3
        if(lmin(i).eq.21) nsphel=nsphel+6
      enddo
      return
      end
c
C*MODULE INT1    *DECK PCMPOT
      SUBROUTINE PCMPOT(H,VPCM,qse,xyzcts,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL IANDJ,DOUBLE,GOPARR,DSKWRK,MASWRK,NXT
      PARAMETER (MXSH=5000, MXGTOT=20000,MXATM=2000,MAXSH=35)
C
      DIMENSION H(L2),VPCM(L2),qse(*),xyzcts(mxts,3),
     *          VPBLK(225),DIJ(225),IJX(225),IJY(225),IJZ(225),
     *          XIN(125),YIN(125),ZIN(125),CONI(MAXSH),CONJ(MAXSH),
     *          IX(35),IY(35),IZ(35),JX(35),JY(35),JZ(35),HP(28),WP(28),
     *          MINP(7),MAXP(7)
C
      COMMON /HERMIT/ H1,H2(2),H3(3),H4(4),H5(5),H6(6),H7(7),
     *                H8(8),H9(9),H10(10)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMDIM/ mxsp,mxts,mempcm1,mempcm2,NTS
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /SHLNRM/ PNRM(84)
c     COMMON /STV   / XINT,YINT,ZINT,TAA,X0,Y0,Z0,
c    *                XI,YI,ZI,XJ,YJ,ZJ,NI,NJ
      COMMON /SYMIND/ TOL,II,JJ,LIT,LJT,MINI,MINJ,MAXI,MAXJ,IANDJ
      COMMON /WERMIT/ W1,W2(2),W3(3),W4(4),W5(5),W6(6),W7(7),
     *                W8(8),W9(9),W10(10)
c
      EQUIVALENCE (HP(1),H1),(WP(1),W1)
C
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00, ONE=1.0D+00, TWO=2.0D+00,
     *           PI212=1.1283791670955D+00, RLN10=2.30258D+00)
C
      DATA MINP /1,2,4,7,11,16,22/
      DATA MAXP /1,3,6,10,15,21,28/
      DATA JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     *          3, 0, 0, 2, 2, 1, 0, 1, 0, 1,
     *          4, 0, 0, 3, 3, 1, 0, 1, 0, 2,
     *          2, 0, 2, 1, 1/
      DATA IX / 1, 6, 1, 1,11, 1, 1, 6, 6, 1,
     *         16, 1, 1,11,11, 6, 1, 6, 1, 6,
     *         21, 1, 1,16,16, 6, 1, 6, 1,11,
     *         11, 1,11, 6, 6/
      DATA JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     *          0, 3, 0, 1, 0, 2, 2, 0, 1, 1,
     *          0, 4, 0, 1, 0, 3, 3, 0, 1, 2,
     *          0, 2, 1, 2, 1/
      DATA IY / 1, 1, 6, 1, 1,11, 1, 6, 1, 6,
     *          1,16, 1, 6, 1,11,11, 1, 6, 6,
     *          1,21, 1, 6, 1,16,16, 1, 6,11,
     *          1,11, 6,11, 6/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     *          0, 0, 3, 0, 1, 0, 1, 2, 2, 1,
     *          0, 0, 4, 0, 1, 0, 1, 3, 3, 0,
     *          2, 2, 1, 1, 2/
      DATA IZ / 1, 1, 1, 6, 1, 1,11, 1, 6, 6,
     *          1, 1,16, 1, 6, 1, 6,11,11, 6,
     *          1, 1,21, 1, 6, 1, 6,16,16, 1,
     *         11,11, 6, 6,11/
C
C     ----- COMPUTE CONVENTIONAL H, S, AND T INTEGRALS -----
C
      call timit(1)
c
      itolpcm=itol-4
      TOL = RLN10*ITOLpcm
      CALL VCLR(vpcm,1,L2)
C
C     ----- INTIALIZE PARALLEL -----
C
      NXT = IBTYP.EQ.1
      NEXT  = -1
      kount = -1
C
C     ----- I SHELL -----
C
      DO 720 II = 1,NSHELL 
         I = KATOM(II)
         XI = C(1,I)
         YI = C(2,I)
         ZI = C(3,I)
         I1 = KSTART(II)
         I2 = I1+KNG(II)-1
         LIT = KTYPE(II)
         MINI = KMIN(II)
         MAXI = KMAX(II)
         LOCI = KLOC(II)-MINI
C
C     ----- J SHELL -----
C
         DO 700 JJ = 1,II
C
C     ----- GO PARALLEL! -----
C
           if(goparr) then 
              kount=kount+1
              if(nxt) then
                 IF(kount.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
                 if(NEXT.ne.kount) goto 700
              else  
                 if(MOD(kount,NPROC).NE.me) goto 700
              endif
            endif
            J = KATOM(JJ)
            XJ = C(1,J)
            YJ = C(2,J)
            ZJ = C(3,J)
            J1 = KSTART(JJ)
            J2 = J1+KNG(JJ)-1
            LJT = KTYPE(JJ)
            MINJ = KMIN(JJ)
            MAXJ = KMAX(JJ)
            LOCJ = KLOC(JJ)-MINJ
            NROOTS = (LIT+LJT-2)/2+1
            RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
            IANDJ = II .EQ. JJ
C
C     ----- PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS
C
            IJ = 0
            JMAX = MAXJ
            DO 160 I = MINI,MAXI
               NX = IX(I)
               NY = IY(I)
               NZ = IZ(I)
               IF (IANDJ) JMAX = I
               DO 140 J = MINJ,JMAX
                  IJ = IJ+1
                  IJX(IJ) = NX+JX(J)
                  IJY(IJ) = NY+JY(J)
                  IJZ(IJ) = NZ+JZ(J)
  140          CONTINUE
  160       CONTINUE
C
            CALL VCLR(VPBLK,1,IJ)
C
C     ----- I PRIMITIVE
C
            JGMAX = J2
            DO 520 IG = I1,I2
               AI = EX(IG)
               ARRI = AI*RR
               CALL SETCONI(CONI,IG)
C
C     ----- J PRIMITIVE
C
               IF (IANDJ) JGMAX = IG
               DO 500 JG = J1,JGMAX
                  AJ = EX(JG)
                  AA = AI+AJ
                  AA1 = ONE/AA
                  DUM = AJ*ARRI*AA1
                  IF (DUM .GT. TOL) GO TO 500
                  FAC = EXP(-DUM)
                  CALL SETCONI(CONJ,JG)
                  AX = (AI*XI+AJ*XJ)*AA1
                  AY = (AI*YI+AJ*YJ)*AA1
                  AZ = (AI*ZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
C     THE ONLY REASON WHY -ILZ WORKS WITH THIS DENSITY THAT ASSUMES
C     HERMITICITY IS BECAUSE <I|-ILZ|I>=0 (MOMENTUM QUENCHING).
C
                  DOUBLE=IANDJ.AND.IG.NE.JG
                  JMAX = MAXJ
                  NN = 0
                  DTWO=ONE
                  IF(DOUBLE) DTWO=TWO
C                 NOTE THAT PNRM FACTORS FOR S AND P SHELLS ARE 1.
                  DUM = PI212*AA1*FAC
                  SPDIJ=CS(IG)*CP(JG)*DUM
                  DO I = MINI,MAXI
                     IF (IANDJ) JMAX = I
                     FACI=DUM*CONI(I)*PNRM(I)*DTWO
                     NN1=NN+1
                     DO J = MINJ,JMAX
                        NN = NN+1
                        DIJ(NN)=FACI*CONJ(J)*PNRM(J)
                     enddo
C            CORRECT FOR L-SHELL DOUBLE COUNTING OF THE SP
C            OFF-DIAGONAL TERMS (FOR NON-L SHELLS CSI*CPJ IS ZERO).
C            NN1 POINTS TO THE APPROPRIATE DENSITY ELEMENT
                     IF(MINJ.LE.1.AND.I.GT.1.AND.DOUBLE)
     *                 DIJ(NN1)=DIJ(NN1)*PT5+SPDIJ
                  enddo 
C
                  AAX = AA*AX
                  AAY = AA*AY
                  AAZ = AA*AZ
C
C     -NCHMAT- IS NONZERO IF THERE ARE EXTERNAL CHARGES WHICH
C     PERTURB THE SYSTEM, SUCH AS IF CHARMM IS IN USE.  NOTE
C     THAT THERE IS ALSO A NUCLEAR REPULSION TERM WHICH IS NOT
C     INCLUDED HERE, IT IS IN THE CHARMM INTERFACE CODE.
C
                  DO 460 IC = 1,NTS
                     ZNUC = -QSE(IC)
                     CX = xyzCTS(IC,1)
                     CY = xyzCTS(IC,2)
                     CZ = xyzCTS(IC,3)
                     if(znuc.eq.zero) goto 460 
c                    if(testtol.and.abs(znuc*fac).lt.thrpcm) goto 460
                     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
                     MM = 0
                     DO 430 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
                        IN = -5+MM
                        DO 320  I=1,LIT
                           IN = IN+5
                           NI = I
C
                           DO 320  J=1,LJT
                              JN = IN+J
                              NJ = J
C
C       EVALUATE MOMENT INTEGRALS USING GAUSS-HERMITE QUADRATURE:
C
                              XINT0 = ZERO
                              YINT0 = ZERO
                              ZINT0 = ZERO
C
                              NPTS = (NI + NJ - 2)/2 + 1
                              IMIN = MINP(NPTS)
                              IMAX = MAXP(NPTS)
C
                              DO 310  IROOT=IMIN,IMAX
C
                                 DUM = WP(IROOT)
                                 PX = DUM
                                 PY = DUM
                                 PZ = DUM
C
                                 DUM = HP(IROOT)*T
                                 PTX = DUM + X0
                                 PTY = DUM + Y0
                                 PTZ = DUM + Z0
C
                                 AXI = PTX - XI
                                 AYI = PTY - YI
                                 AZI = PTZ - ZI
C
                                 BXI = PTX - XJ
                                 BYI = PTY - YJ
                                 BZI = PTZ - ZJ
C
                                 GO TO (250,240,230,220,210),NI
C
  210                            PX = PX*AXI
                                 PY = PY*AYI
                                 PZ = PZ*AZI
C
  220                            PX = PX*AXI
                                 PY = PY*AYI
                                 PZ = PZ*AZI
C
  230                            PX = PX*AXI
                                 PY = PY*AYI
                                 PZ = PZ*AZI
C
  240                            PX = PX*AXI
                                 PY = PY*AYI
                                 PZ = PZ*AZI
C
  250                            CONTINUE
C
                                 GO TO (300,290,280,270,260),NJ
C
  260                            PX = PX*BXI
                                 PY = PY*BYI
                                 PZ = PZ*BZI
C
  270                            PX = PX*BXI
                                 PY = PY*BYI
                                 PZ = PZ*BZI
C
  280                            PX = PX*BXI
                                 PY = PY*BYI
                                 PZ = PZ*BZI
C
  290                            PX = PX*BXI
                                 PY = PY*BYI
                                 PZ = PZ*BZI
C
  300                            CONTINUE
C
                                 XINT0 = XINT0 + PX
                                 YINT0 = YINT0 + PY
                                 ZINT0 = ZINT0 + PZ
C
  310                         CONTINUE
C
                              XIN(JN) = XINT0
                              YIN(JN) = YINT0
                              ZIN(JN) = ZINT0*WW
  320                   CONTINUE
                        MM = MM+25
  430                CONTINUE
C
                     DO 450 I = 1,IJ
                        NX = IJX(I)
                        NY = IJY(I)
                        NZ = IJZ(I)
                        DUM = ZERO
                        MM = 0
                        DO 440 K = 1,NROOTS
                           DUM = DUM+XIN(NX+MM)*YIN(NY+MM)*ZIN(NZ+MM)
                           MM = MM+25
  440                   CONTINUE
                        VPBLK(I) = VPBLK(I) + DUM*DIJ(I) 
  450                CONTINUE
C
  460             CONTINUE
C
C     ----- END OF PRIMITIVE LOOPS -----
C
  500          CONTINUE
  520       CONTINUE
C
C     ----- COPY BLOCK INTO H-CORE, OVERLAP, AND KINETIC ENERGY MATRICES
C
            JMAX = MAXJ
            NN = 0
            DO 620 I = MINI,MAXI
               LI = LOCI+I
               IN = (LI*(LI-1))/2
               IF (IANDJ) JMAX = I
               DO 600 J = MINJ,JMAX
                  LJ = LOCJ+J
                  JN = LJ+IN
                  NN = NN+1
                  VPCM(JN)=VPBLK(NN)
  600          CONTINUE
  620       CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
  700    CONTINUE
  720 CONTINUE
C
C     ----- SUM UP PARTIAL CONTRIBUTIONS IF PARALLEL -----
C
      IF (GOPARR) then
         IF(nxt) CALL DDI_DLBRESET
         CALL DDI_GSUMF(914,VPCM,L2)
      endif
c
      CALL DAWRIT(IDAF,IODA,VPCM,L2,311,0)
c     ADD VPCM to the core (aka 1e) Hamiltonian.
      CALL DAread(IDAF,IODA,H,L2,11,0)
      call daxpy(l2,one,VPCM,1,H,1)
      CALL DAWRIT(IDAF,IODA,H,L2,11,0)
      if(maswrk) write(iw,*) 'Done VPCM.'
      call timit(1)
c     write(6,*) 'wwwpcm-int'
c     call prtrila(VPCM,num)
      RETURN
      END
C
