C 22 Feb 06 - TN  - MAKMOL: read $EFRAG for EFP/FMO model
C 14 Nov 05 - DGF - various changes for FMO 2.1 release
C 19 Sep 05 - MWS - add true nuclear charge array to INFOA common
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 - mws - pad common FMORUN and common MCINP
C 23 Jul 04 - MWS - FMOINP: frgnam is D.P., not character type
c 19 May 04 - DGF,KK - implement Fragment Molecular Orbital (FMO) method
c
C*MODULE fmolib  *DECK makemol
      SUBROUTINE makemol(ifg,jfg,kfg,ilayer,itask,nat0,nshell0,ngau0,ne0
     *                  ,ich0,mul0,runqm)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical runqm
      COMMON /FMCOM / X(1)
      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
c
      call makmol(ifg,jfg,kfg,ilayer,itask,nat0,nshell0,ngau0,ne0,ich0,
     *            mul0,runqm,x(lichfg),x(lmulfg),x(lfrgnam),x(liaglob),
     *            x(liatfrg),x(lindfrg),x(lindgfrg),x(lnatfrg),
     *            x(lnat0frg),x(lianfrg),x(lzanfrg),x(lcfrg),x(llibish),
     *            x(llibnsh),x(llibng),x(lizbas),x(lloctat),x(liaoglob))
c
      return
      end
C*MODULE fmolib  *DECK makmol
      SUBROUTINE makmol(ifg,jfg,kfg,ilayer,itask,nat0,nshell0,ngau0,ne0,
     *                  ich0,mul0,runqm,ichfg,mulfg,frgnam,iaglob,iatfrg
     *                ,indfrg,indgfrg,natfrg,nat0frg,ianfrg,zanfrg,cfrg,
     *                  libish,libnsh,libng,izbas,loctat,iaoglob)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical GOPARR,DSKWRK,MASWRK,isgddi,parout,some,mptest,runqm,
     *        totprop,LEFNUC
c
      CHARACTER*4  ATM
      character*8 ATOMNM,frgnam(*),c8dum
      character*4 atomnm4
      Parameter (MXATM=2000,MXSH=5000,MXGTOT=20000,MXAO=8192,MXRT=100,
     *           MAXNZ=137,MAXL=5,MaxLay=5,
     *           MXSFMO=MaxLay*40,MXGFMO=MaxLay*100,MXAFMO=MaxLay*10)
      dimension ichfg(*),mulfg(*),iaglob(*),iatfrg(*),indfrg(*),
     *          indgfrg(*),natfrg(*),nat0frg(*),ianfrg(*),zanfrg(*),
     *          cfrg(3,*),libish(MAXNZ,maxbas,*),libnsh(MAXNZ,maxbas,*),
     *          libng(MAXNZ,maxbas,*),izbas(*),loctat(*),iaoglob(*)
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /CORE  / CORE(107)
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /EFPFMO/ IEFFMO,LEFNUC,REPNUCEFP
      COMMON /GDDI/   ISCOPE,NGROUPS,MYGROUP,MEGLOB,npglob,nnglob,JBTYP,
     *                ISGDDI,PAROUT
      COMMON /INFO  / CDUM(MXATM,3),IXXX(MXATM),NATOMS,IUNTRD,ATM(MXATM)
      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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SOOPT / NOSO
      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 /fmorun/ espscf,IDAFMO,icurlay,icurfg,jcurfg,kcurfg,ncursh,
     *                ngau,icurpop,ifmostp,moncor,nat1e,needr,modrst,
     *                norbproj
      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
      equivalence (ATOMNM,FATOMNM)
      data blank/8H        /,none/4HNONE/,rnone/8HNONE    /,
     *     RMC/8HMCSCF   /,dbgfmo/8HDBGFMO  /,dbgme/8HMAKMOL  /,
     *     debug/8HDEBUG   /,RHF/8HRHF     /,CIS/8HCIS     /,c8dum/' '/
c
c     itask=0 make coordinates+basis set
c     itask=1 make coordinates (light version) 
c     parstat: GroupNone
c
      some=(exetyp.eq.debug.or.exetyp.eq.dbgfmo.or.exetyp.eq.dbgme).and.
     *     maswrk
      totprop=iand(modprp,1).ne.0.and.runqm
c
c     a) create atomic coordinates
c
      nat=nat0
      nati=nat0frg(ifg)
      if(jfg.eq.0) nati=natfrg(ifg)
      indi=indfrg(ifg)
      call addfrg(nati,indi,zanfrg,cfrg,ianfrg,iatfrg,iaglob)
      nati=natfrg(ifg)
      nat0i=nat0frg(ifg)
      indgi=indgfrg(ifg)
c     first add nonghost atoms. Note that if only monomer is requested we add
c     also ghost atoms since they are never redundant.
      natj=0
      nat0j=0
      indj=1
      indgj=1
      natk=0
      nat0k=0
      indk=1
      indgk=1
      if(jfg.ne.0) then
        indj=indfrg(jfg)
        natj=natfrg(jfg)
        nat0j=nat0frg(jfg)
        indgj=indgfrg(jfg)
        call addfrg(nat0j,indj,zanfrg,cfrg,ianfrg,iatfrg,iaglob)
        if(kfg.ne.0) then
          indk=indfrg(kfg)
          natk=natfrg(kfg)
          nat0k=nat0frg(kfg)
          indgk=indgfrg(kfg)
          call addfrg(nat0k,indk,zanfrg,cfrg,ianfrg,iatfrg,iaglob)
        endif
      endif
c     add ghost atoms. Some can be redundant with real atoms so one has to
c     sift through (redundant if a broken bond is inside the n-mer).
      if(nbdfg.ne.0) then 
        if(jfg.ne.0) then
        natgi=nati-nat0i
        natgj=natj-nat0j
        natgk=natk-nat0k
c       sift ghost atoms of ifg in jfg and kfg 
        if(natgi.ne.0) call addgho(natgi,zanfrg(indgi),cfrg(1,indgi),
     *     ianfrg(indgi),iatfrg(indgi),nat0j,iatfrg(indj),nat0i,nat0k,
     *     iatfrg(indk),nat0i+nat0j,iaglob)
c       sift ghost atoms of jfg in ifg and kfg 
        if(natgj.ne.0) call addgho(natgj,zanfrg(indgj),cfrg(1,indgj),
     *     ianfrg(indgj),iatfrg(indgj),nat0i,iatfrg(indi),0,    nat0k,
     *     iatfrg(indk),nat0i+nat0j,iaglob)
c       sift ghost atoms of kfg in ifg and jfg 
        if(natgk.ne.0) call addgho(natgk,zanfrg(indgk),cfrg(1,indgk),
     *     ianfrg(indgk),iatfrg(indgk),nat0i,iatfrg(indi),0,    nat0j,
     *     iatfrg(indj),nat0i,      iaglob)
        endif
c       do ibdfg=1,nbdfg
c         ibdgh(ibdfg)=0
c         iafg=indat(abs(iabdfg(ibdfg)))
c         jafg=indat(jabdfg(ibdfg))
c         if(iafg.ne.ifg.and.iafg.ne.jfg.and.iafg.ne.kfg.and.
c    *       (jafg.eq.ifg.or.jafg.eq.jfg.or.jafg.eq.kfg)) 
c    *      ibdgh(ibdfg)=iand(ialoc(natfmo+ibdfg),65535)
c       enddo
      endif
      if(some) write(6,6666) (zan(kk),(c(jj,kk),jj=1,3),kk=1,nat)
 6666 format(100(F10.5,3F12.5,/))
c
c     b) copy basis set from the library
c
      ne=ne0
      ii=nshell0
      ish=ngau0
c     write(6,*) 'wwwsh',nshell0,ngau0
      do iat=nat0+1,nat
        izi=ian(iat)
c       get the basis set for this atom. Use the atomic number, not charge.
        ibas=izbas(iaglob(iat))
        jj=libish(izi,ibas,ilayer)
        if(jj.eq.0) then
          write(iw,9120) izi,ibas,iat,ifg,ilayer
          call abrt
        endif
        nsh=libnsh(izi,ibas,ilayer)
        ngg=libng(izi,ibas,ilayer)
c       write(6,*) 'wwwshh',nsh,ngg,jj,izi,ibas,ilayer
        if(ii+nsh.gt.MXSH) then
          if(maswrk) write(iw,*) 'Increase MXSH',ii+nsh,MXSH
          call abrt
        endif
        if(ish+ngg.gt.MXGTOT) then
          if(maswrk) write(iw,*) 'Increase MXGTOT',ish+ngg,MXGTOT
          call abrt
        endif
        jsh=lstart(jj)
        if(itask.eq.0) then
          do k=1,ngg
            ish=ish+1
            EX(ish)=fex(jsh) 
            CS(ish)=fC(jsh,1)
            CP(ish)=fC(jsh,2)
            CD(ish)=fC(jsh,3)
            CF(ish)=fC(jsh,4)
            CG(ish)=fC(jsh,5)
            jsh=jsh+1
          enddo
        endif
        if(totprop) iaoff=loctat(iaglob(iat))
        do i=1,nsh
          ii=ii+1
          kATOM(ii) =iat
          kTYPE(ii) =lTYPE(jj)
          kNG(ii)   =lNG(jj)
          if(ii.eq.1) then
            kstart(ii)=1
            kloc(ii)=1
          else
            kSTART(ii)=kstart(ii-1)+kng(ii-1)
            kloc(ii)=kloc(ii-1)+(kmax(ii-1)-kmin(ii-1)+1)
          endif
          kMIN(ii)  =lMIN(jj)
          kMAX(ii)  =lMAX(jj)
          jj=jj+1
          if(totprop) then
c           This code for building AOi map is not very safe if basis set
c           storage changes. In particular, it is assumed here that all 
c           shells for each atom are grouped together.
            iaomin=kloc(ii)
            iaomax=iaomin+kMAX(ii)-kMIN(ii)
            do iao=iaomin,iaomax
              iaoglob(iao)=iaoff
              iaoff=iaoff+1
            enddo
          endif
        enddo
        atomnm4='X '
        imode=ian(iat)
        if(imode.ne.0) call zsymnum(c8dum,atomnm4,imode)
        write(UNIT=ATOMNM,FMT='(A2,I6)') atomnm4,iat
        anam(iat)=FATOMNM
        bnam(iat)=BLANK
c       what is BNAM??
c       nei=ian(iat)
        izan=int(zan(iat)+0.5D+00)
        iian=ian(iat)
        nei=izan
c       compensate for double counting of charges of broken bonds.
        nei=izan
        if(izan.ne.iian) then
          if(izan.eq.1) then
c           nei=nei-1 
          else
c           nei=nei+1 
          endif
        endif
        if(MPCTYP.ne.NONE.and.imode.ne.0) nei=INT(CORE(imode))
        ne=ne+nei
      enddo
      nshell=ii
      num=kloc(ii)+kmax(ii)-kmin(ii)
      nqmt=num
      ich=ich0+ichfg(ifg)
      mul=1
      if(scftyp.ne.rhf) then
        mul=mulfg(ifg)
c       change multiplicity only for non-RHF (at present really only MCSCF). 
c       in case of RHF one can ask for CI of other multiplicities while
c       getting the orbitals from RHF for singlets. In this case CI
c       multiplicities will be taken from mulfg.
        if(mul0.ne.0) then
c          Allow at most one non-singlet
           if(mul0.ne.1.and.mul.ne.1) then
             if(maswrk) write(iw,9110) 'ESP',mul,mul0,ifg,jfg
             call abrt
           endif
           mul=max(mul,mul0)
        endif
      endif
      if(jfg.ne.0) then
        ich=ich+ichfg(jfg)
        if(scftyp.ne.rhf) then
          mulj=mulfg(jfg)
c         Allow at most one non-singlet
          if(mul.ne.1.and.mulj.ne.1) then
            if(maswrk) write(iw,9110) 'dim',mul,mul0,ifg,jfg
            call abrt
          endif
          mul=max(mul,mulj)
        endif
      endif
      if(kfg.ne.0) then
        ich=ich+ichfg(kfg)
        if(scftyp.ne.rhf) then
          mulk=mulfg(kfg)
c         Allow at most one non-singlet
          if(mul.ne.1.and.mulk.ne.1) then
            if(maswrk) write(iw,9110) 'trim',mul,mul0,ifg,jfg,kfg
            call abrt
          endif
          mul=max(mul,mulk)
        endif
      endif
      ngau=ish
c
      NE = NE-ICH
      NA=(NE+MUL-1)/2
      NB=(NE-MUL+1)/2
c
c     The check below is stolen from ATOMS.
c 
      IF(NA+NB .NE. NE) THEN
         IF (MASWRK) WRITE(IW,9290) ifg,NE,ICH,MUL
         CALL ABRT
      END IF
C           IMPOSSIBLY HIGH, LOW, OR MISMATCHED TO E- COUNT
      IDUM=MUL+NE
      IF(MUL.GT.NE+1 .OR. MUL.LT.0 .OR. 2*INT(IDUM/2).EQ.IDUM) THEN
        IF (MASWRK) WRITE(IW,9280) ifg,mul,ne
        CALL ABRT
      END IF
c
      MPCSAV = MPCTYP
      MPCTYP = NONE
      ENUCR  = ENUC(NAT,ZAN,C)
      MPCTYP = MPCSAV
      CALL SETLAB(1,atomnm4)
cnb5  This should be adjusted for MCQDPT.
c     MP2inp sets up some variables so we have to call it with every new
c     fragment. ONly call MP2inp when external monomers are not around.
      MPTEST=.FALSE.
      if(itask.eq.0.and.(icurfg.eq.0.or.ifmostp.eq.1.and.nbsse.eq.3)
     *   .and.ifmostp.ne.6) then
        if(mplevl.ne.0) call MP2INP(mptest)
        IF(CCTYP.NE.RNONE) CALL CCINP
        IF(SCFTYP.EQ.RMC) CALL MCIN
        IF(CITYP.EQ.CIS) CALL CISINP
        IF(IEFFMO.EQ.1) CALL EFINP(IUNTRD,IEF)
      endif
c     if(itask.eq.0.and..not.dirscf.and..not.dirtrf) call mod2ei
      if(itask.eq.0) call mod2ei
c
c     finally, fill in symmetry (force C1).
c
c     nt=1
c     t(1)=one
c     invt(1)=1
c
c     ispher=1 means QMTSYM must eliminate extra functions. 
c     With RHF we try to save time by not getting the Q-matrix but
c     MCSCF apparently needs one even for C1.
c
      if(itask.eq.0.and.some) then
         if(jfg.eq.0) then
           write(iw,9010) frgnam(ifg),ifg,mygroup
         else if(kfg.eq.0) then
           write(iw,9020) frgnam(ifg),frgnam(jfg),ifg,jfg,mygroup
         else
           write(iw,9030) frgnam(ifg),frgnam(jfg),frgnam(kfg),ifg,jfg,
     *                    kfg,mygroup
         endif
         write(iw,9000) nshell,num,ne,ich,mul,na,nb,nat,enucr
         call intr
      endif
c     icoord is not known
c     IF(runqm.and.ICOORD.NE.4) THEN
      IF(runqm) THEN
c        jrest=irest
         if(irest.gt.1) irest=0
c        if SCF does not converge IREST is set to 2. This prevents further
c        monomers/dimers from running by not saving SPDTR matrices
c        SPDTR matrices are used to rotate orbital coefficients 
         if(nbdfg.ne.0.or.ispher.ge.0.or.SCFTYP.eq.rmc) CALL SPDTR
         if(ispher.lt.0) noso=1
         if(scftyp.eq.rmc) noso=0
         if(ispher.ge.0.or.scftyp.eq.rmc) CALL SYMORB
      END IF
c     if(totprop) write(6,*) 'wwwiaglob',(iaoglob(i),i=1,num)
      return
 9010 format(/1x,'Monomer fragment ',A8,'(',I3,') is done by the group',
     *           I3/)
 9020 format(/1x,'Dimer fragment ',2A8,'(',2I3,') is done by the group',
     *           I3/)
 9030 format(/1x,'Trimer fragment ',3A8,'(',3I3,') is done by the group'
     *          ,I3/)
 9000 format(/1X,'TOTAL NUMBER OF BASIS SET SHELLS             =',I5/
     *        1X,'NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =',I5/
     *        1X,'NUMBER OF ELECTRONS                          =',I5/
     *        1X,'CHARGE OF MOLECULE                           =',I5/
     *        1X,'SPIN MULTIPLICITY                            =',I5/
     *        1X,'NUMBER OF OCCUPIED ORBITALS (ALPHA)          =',I5/
     *        1X,'NUMBER OF OCCUPIED ORBITALS (BETA )          =',I5/
     *        1X,'TOTAL NUMBER OF ATOMS                        =',I5/
     *        1x,'THE NUCLEAR REPULSION ENERGY IS',F20.10)
 9110 format(/1x,'Check multiplicities, at most one nonsinglet allowed:'
     ,          ,A3,1x,2I2,2I4/)
 9120 format(/1x,'Basis set not found for iz=',I4,', ibas=',i3,', iat=',
     *           I6,', ifg=',I5,', ilay=',i2,'.',/)
 9280 FORMAT(/1X,'Impossible spin multiplicity',I3,' with',I6,
     *           ' electrons for fragment',I6,/)
 9290 FORMAT(/1X,'Fragment',I6,':',I6,' electrons, charge',I3,
     *           ', multiplicity',I3,' - impossible!',/)
      end
c
C*MODULE fmolib  *DECK makesbs
      SUBROUTINE makesbs(ij,ifg,zsave,nesav,indat,iaglob)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      Parameter (MXATM=2000,zero=0.0D+00)
      logical mptest
      dimension zsave(*),indat(*),iaglob(*)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      data rnone/8HNONE    /,RMC/8HMCSCF   /,CIS/8HCIS     /
c
c     prepare basis sets and other info for BSSE corrections.
c     zero out charges for BSSE ghost atoms.
c     this subroutines is called twice, 
c     ij=1 : zero out charges for fragment J (must be called before ij=2)
c     ij=2 : zero out charges for fragment I 
c
      if(ij.eq.1) then
         call dcopy(nat,zan,1,zsave,1)
         nesav=ne
      else
         call dcopy(nat,zsave,1,zan,1)
         ne=nesav
      endif
c
cnb5  adjust border charges
cnb5  adjust multiplicity 
c
      nei=0
      do i=1,nat
        if(indat(iaglob(i)).ne.ifg) then
           nei=nei+int(zan(i)+0.5D+00)
           zan(i)=zero
        endif
      enddo
      NE = ne - nei
      NA = (NE+MUL-1)/2
      NB = (NE-MUL+1)/2
c
c     initialise correlation runs (under construction!). 
c
      MPTEST=.FALSE.
      if(mplevl.ne.0) call MP2INP(mptest)
      IF(CCTYP.NE.RNONE) CALL CCINP
      IF(SCFTYP.EQ.RMC) CALL MCIN
      IF(CITYP.EQ.CIS) CALL CISINP
c
      RETURN
      END
C*MODULE fmolib  *DECK projgues
      SUBROUTINE projgues(l1,l2,vv,wrk1,wrk2,wrk3)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      dimension vv(*),wrk1(l1),wrk2(l1*l1),wrk3(l2) 
      PARAMETER (MXSH=5000,MXGTOT=20000)
      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 /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /FMCOM / X(1)
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      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
c
c
c     Project out cut bond LMOs. This is useful in making initial guess for
c     dimers using monomer orbitals.
c
c     do nothing if not exchanging orbitals
c
      if(mod(modorb,2).eq.0) return
c
c     l3=l1*l1
c     CALL DAREAD(IDAF,IODA,vv,L3,15,0)
      CALL VALFM(LOADFM)
      lss=LOADFM+1
      ldd=lss+l1*l1
      lrotlcao=ldd+(maxcbs*maxcbs+maxcbs)/2
      last=lrotlcao+maxcbs*maxcao
      lq=lss
      lscr=lss
      NEED = LAST- LOADFM -1
      CALL GETFM(NEED)
      CALL DAread(IDAF,IODA,wrk3,L2,12,0)
      call fmohop(l1,l2,vv,wrk3,x(lss),x(lq),x(ldd),x(lscr),wrk2,wrk1,
     *            x(liabdfg),x(ljabdfg),x(lidxCAO),x(liaglob),x(lnCBS),
     *            x(lnCAO),x(liaprjo),x(ljaprjo),x(lshiftb),x(lCoreAO),
     *            x(lfmoc),x(lrotlcao),x(llocfmo),nshell,KATOM,KTYPE,
     *            KLOC,kmin,.false.)
c     call prsq(vv,l1,l1,l1)
c     call TFTRI(wrk2,da,wrk1,x(lwrk),l1,l1,l1)
c     call dcopy(l2,wrk1,1,da,1)
      CALL RETFM(NEED)
c     CALL dawrit(IDAF,IODA,vv,L3,15,0)
      if(maswrk) write(iw,*) 'Orbitals have been projected'
      RETURN
      END
C*MODULE fmolib  *DECK madtrap
      SUBROUTINE madtrap(m,a,n,ix,b,c)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      dimension a(*),b(*),c(n,n)
c
c     C=C+A*B, where A is triangular, B is trapezoid (part of a triangular
c     matrix with linear dimension n starting at location ix) and C is square.
c     (only a block of C that has size m*n is updated by addition)
c
      if(m.eq.0.or.n.eq.0.or.m+ix-1.gt.n) then
        write(6,*) 'bad indices',m,n,m+ix
        return
      endif
      do i=1,m
        ic=ix+i-1
        do j=1,n
          sum=0.0D+00
          do k=1,m
            ia=max(i,k)
            ia2=(ia*ia-ia)/2
            ib=ix+k-1
            ib2=max(ib,j)
            sum=sum+a(ia2+min(i,k))*b((ib2*ib2-ib2)/2+min(ib,j))
          enddo
          c(ic,j)=c(ic,j)+sum
        enddo
      enddo
      RETURN
      END
C*MODULE fmolib  *DECK fmoatfrg
      SUBROUTINE fmoatfrg(iat0,indat,indatg,iaglob,ialoc,iabdfg,jabdfg,
     *                    indbd,fmozan,fmoc,popmat,iats,jats,fracv,
     *                    fracesp,iatfrg,ZNUC,cx,cy,cz)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      logical iainside,jainside,smartr(2),bimer(3)
      dimension indat(*),indatg(natfmo,*),iaglob(*),iabdfg(*),jabdfg(*),
     *          indbd(maxabd,*),fmozan(*),fmoc(3,*),popmat(maxnat,nfg,2)
     *         ,ialoc(*)
      Parameter (MXATM=2000, zero=0.0D+00,one=1.0D+00,half=0.5D+00)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      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 bimer/3*.false./
c
c     Determine if atom IAT belongs to fragment 
c
      if(ifmostp.ne.6) then
        iat=iat0
        znuc=-fmozan(iat)
        cx=fmoc(1,iat)
        cy=fmoc(2,iat)
        cz=fmoc(3,iat)
        iz=1
        ifg=icurfg
        factk=one
      else
        iat=iat0+nat
        znuc=-zan(iat)
        cx=c(1,iat)
        cy=c(2,iat)
        cz=c(3,iat)
c       iat below must be global iat
        iat=iaglob(iat)
        iz=2
        ifg=ncursh
        factk=half
      endif  
      jfg=jcurfg
      lfg=kcurfg
      iatfrg=indat(iat)
      smartr(1)=iand(ixesp,8).ne.0.and.jfg.ne.0
      smartr(2)=iand(ixesp,16).ne.0.and.jfg.ne.0
      if(smartr(1).and.nbdfg.ne.0) then
        bimer(1)=fmodist(ifg,0,0,jfg).eq.0
        if(lfg.eq.0) then
          if(bimer(1)) smartr(1)=.false.
        else
          bimer(2)=fmodist(ifg,0,0,lfg).eq.0
          bimer(3)=fmodist(jfg,0,0,lfg).eq.0
        if(bimer(1).and.(bimer(2).or.bimer(3)).or.bimer(2).and.bimer(3))
     *    smartr(1)=.false.
        endif 
      endif 
      iatsg=iaglob(iats)
      jatsg=iaglob(jats)
      iifg=indat(iatsg)
      jjfg=indat(jatsg)
c     3 -body terms do not work with BSSE.
      if(nbdfg.eq.0.or.ifmostp.eq.6) then
        fracv=zero
        fracesp=one
c       BSSE
        if(ifmostp.eq.5) then
          if(nbsse.eq.1.and.iatfrg.eq.jfg .or. nbsse.eq.2.and.
     *      iatfrg.ne.ifg.and.(iifg.eq.jfg.or.jjfg.eq.jfg)) fracesp=zero
        endif
c       decide if the point charge iat should be included into ESP (keep=1)
        if(ifmostp.ne.6) then
         if(fracesp.eq.one.and.(nbsse.ne.2.or.ifmostp.ne.5.or.iatfrg.ne.
     *  jfg).and.(iatfrg.eq.ifg.or.iatfrg.eq.jfg.or.iatfrg.eq.lfg)) then
          fracv=one
          fracesp=zero
         endif 
        endif
      else
c       this code does not support nbsse=1,2 and 3.
        if(nbsse.ne.0) then
          write(6,*) 'bad nbsse',nbsse
          call abrt
        endif
c       set values for the case atom IAT is not shared
        if(iatfrg.eq.ifg.or.iatfrg.eq.jfg.or.iatfrg.eq.lfg) then
          fracv=one
          fracesp=zero
        else
          fracv=zero
          fracesp=one
        endif
        fracesp1=zero
c       do ibdfg=1,nbdfg
c       loop over all broken bonds in which atom IAT is involved.
c       we are only interested really in one n-mer here (n=1,2), so
c       we find if an atom is split between this n-mer and the rest, then quit. 
        do iabd=1,maxabd
          ibdfg=indbd(iabd,iat)
          if(ibdfg.eq.0) goto 100
          ia=abs(iabdfg(ibdfg))
          ja=abs(jabdfg(ibdfg))
c         reordering is now done in fmobon
c         the code below assumes the canonical order (left is negative). 
c         if(ja.lt.0) then
c           if(ia.lt.0) call abrt
c           idum=ia
c           ia=ja
c           ja=idum
c         endif
c         ia=abs(ia)
          iafrg=indat(ia)
          jafrg=indat(ja)
          iainside=iafrg.eq.ifg.or.iafrg.eq.jfg.or.iafrg.eq.lfg
          jainside=jafrg.eq.ifg.or.jafrg.eq.jfg.or.jafrg.eq.lfg
c         exclude the case when a cut bond is inside a dimer
          if(iat.eq.ia.and..not.(iainside.and.jainside)) then
            if(iainside) then
              fracv=(znuc+one)/znuc
              fracesp=-one/znuc
            endif
            if(jainside) then
              fracv=-one/znuc
              fracesp=(znuc+one)/znuc
            endif
          endif
c
c           now add partial charges. 
c           Here we add the "right" side (that is, pseudoproton). 
c           The left side (N-1) is treated below, along with all other cases. 
c           
          if(iat.eq.ia.and..not.jainside.and.resppc(iz).ne.zero) then
c           kfg is equal to jafrg.
c           tricky part! kat comes from the ghost atom stored in ialoc
c           Note that this atom has coordinates equal to those of iat
c           but its charge is stored as processed below.
            kfg=ishft(ialoc(natfmo+ibdfg),-16)
            kat=iand(ialoc(natfmo+ibdfg),65535)
cnb         ifmostp.eq.6 cannot come here?!
            if(ifmostp.ne.6) then
              if(smartr(1)) then
                rk=fmossdist(iifg,jjfg,indatg(iatsg,1),indatg(jatsg,1),
     *                      ifg,jfg,lfg,kfg,bimer)
              else if(smartr(2)) then
                rk=fmosdist(iifg,jjfg,indatg(iatsg,1),indatg(jatsg,1),
     *                      ifg,jfg,lfg,kfg)
              else
                rk=fmodist(ifg,jfg,lfg,kfg)
              endif
            else
              rk=fmodist(ifg,0,0,kfg)
            endif
            if(rk.gt.resppc(iz)) then 
              fracesp1=fracesp1+popmat(kat,kfg,icurpop)*factk/znuc
            endif 
c           goto 100
          endif
          if(ifmostp.eq.6.and.resppc(iz).ne.zero) then
            write(6,*) 'fmoptc is not programmed for resppc(2) yet'
            call abrt
          endif
        enddo
  100   continue
        fracesp=fracesp+fracesp1
      endif
c     Add approximate 2e ESP contributions that become atomic Mulliken charges
c     centred at the atomic coordinates (point charge approximation).
      if(resppc(iz).ne.zero) then
        ind=iat
        kfg=ishft(ialoc(ind),-16)
        kat=iand(ialoc(ind),65535)
c       check if the gun is accidently or mistakenly loaded.  
        if(kfg.eq.0) call abrt 
c         
c       for regular runs we want all charges except those from I and J. 
        if(ifmostp.ne.6.and.(ifg.eq.kfg.or.jfg.eq.kfg.or.lfg.eq.kfg)) 
     *    goto 200
c       for esdim we only want charges coming from J
c
        if(ifmostp.ne.6) then
          if(smartr(1)) then
            rk=fmossdist(iifg,jjfg,indatg(iatsg,1),indatg(jatsg,1),
     *                   ifg,jfg,lfg,kfg,bimer)
          else if(smartr(2)) then
            rk=fmosdist(iifg,jjfg,indatg(iatsg,1),indatg(jatsg,1),
     *                  ifg,jfg,lfg,kfg)
          else
            rk=fmodist(ifg,jfg,lfg,kfg)
          endif
        else
          rk=fmodist(ifg,0,0,kfg)
        endif
        if(rk.gt.resppc(iz)) then 
          fracesp=fracesp+popmat(kat,kfg,icurpop)*factk/znuc
        endif 
  200   continue
      endif
      RETURN
      END
C*MODULE fmolib  *DECK DMTX2
      SUBROUTINE DMTX2(D,V,M,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION D(*),V(NDIM,M)
C
      PARAMETER (ZERO=0.0D+00,two=2.0D+00)
c
c     this is an RHF version of DMTX. It should/will be replaced later
c     after generalising FMO into ROHF/UHF.
C
      IJ = 0
      DO 120 I = 1,N
         DO 110 J = 1,I
            IJ = IJ + 1
            DUM = ZERO
            DO 100 K = 1,M
               DUM = DUM+V(I,K)*V(J,K)
  100       CONTINUE
            D(IJ) = two*DUM
  110    CONTINUE
  120 CONTINUE
      RETURN
      END
c
C     *MODULE fmolib  *DECK indsort
      subroutine indsort(n,ia,ind)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      integer ia(n),ind(n)
c
c     Find element indices corresponding to decreasing order in array IA.
c     example: input  IA=10,20,15,25,12
c              output IND=4,2,3,5,1
c     IA is not destroyed on exit.
c     This subroutine was provided by T. Ikegami (AIST).
c
      do i = 1, n
         ind(i) = i
      end do

C     Make heap
C      
      do j = (n + 1) / 2, 1, -1
         jj = j
         i = ind(jj)
 10      k = jj * 2
         if (k.le.n) then
            if ( k.ne.n .and. ia(ind(k)).gt.ia(ind(k+1)) ) k = k+1
            if ( ia(i).le.ia(ind(k)) ) goto 20
            ind(jj) = ind(k)
            jj = k
            goto 10
         end if
 20      ind(jj) = i
      end do
      
C     Do sort -- ind(1) is the index to the minimum here.
C
      do m = n-1, 1, -1
         i = ind(m + 1)
         ind(m + 1) = ind(1)
         jj = 1
 30      k = jj * 2
         if (k.le.m) then
            if (k.ne.m .and. ia(ind(k)).gt.ia(ind(k+1))) k = k+1
            if ( ia(i).le.ia(ind(k)) ) goto 40
            ind(jj) = ind(k)
            jj = k
            goto 30
         end if
 40      ind(jj) = i
      end do
      end
C*MODULE fmolib  *DECK tribrk
      subroutine tribrk(ind,is,i,j)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (half=0.5D+00,tiny=1.0D-06)
c
c     break combined triangular index into two linear indices.
c     that is, find such i and j so that
c     ind=i*(i+is)/2+j 
c     input: ind,is (is is either +1 or -1, ind must define positive i,j)
c     output: i'=i+1,j
c     note: for triangular matrices normally is=-1.
c     Nota bene: i is produced raised by one, that is, the subroutine is suited
c     really for is=-1, assuming the diagonal is not stored.
c 
      i=int(-is*half+sqrt(half*is*is+2*ind)+tiny)
      j=ind-(i*i+is*i)/2
      if(j.eq.0) then
         i=i-1
         j=i
      endif
      i=i+1
c     
      return
      end
c
C*MODULE fmolib  *DECK cubbrk
      subroutine cubbrk(ind,is,i,j,k)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (one=1.0D+00,three=3.0D+00,tiny=1.0D-06)
c
c     break combined pyramidal index into three linear indices.
c     that is, find such i,j and k so that
c     ind=(i-1)(i-2)(i-3)/6+j(j+1)/2+k
c     input: ind
c     output: i,j,k
c     note: at present only is=-1 is suported.
c     the equation below is obtained by analytically solving a cubic equation
c     (hint: Mathematica can do it algebraically).
c
c     if(is.ne.-1) call abrt
      a=ind*162.0D+00
      d=((a+sqrt(a*a-108.0D+00))/2.0D+00)**(one/three)
      i=int(one+one/d+d/three+tiny)
      indij=ind-(i*(i-1)*(i-2))/6
      if(indij.ne.0) then
        i=i+1
      else
        indij=indij+((i-1)*(i-2))/2
      endif  
      call tribrk(indij,is,j,k)
      return
      end
c
C*MODULE fmolib  *DECK fillind
      subroutine fillind(nfg2,nfg3,loadhf,numfrg,iwrk,loadm,loadd,loadt,
     *                   maxl30,layfrg,indat,indatg,nCBS,iabdfg,jabdfg,
     *                   idxcao,scffrg,ichfg,ichfmo,nefmo,l0fmo,l1fmo,
     *                   orbxch,enexch,docas,loadbf,loadgr,nstjob,
     *                   semidyn,spargrid,noffg,indgrd,vdwrad,grdpad,
     *                   nxyzg,maxg,some)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      Parameter (MXATM=2000)
      logical orbxch,enexch,docas,semidyn,spargrid,some
      dimension numfrg(*),iwrk(*),loadm(*),loadd(*),loadt(*),maxl30(*),
     *          layfrg(*),indat(*),indatg(natfmo,*),iabdfg(*),jabdfg(*),
     *          idxcao(MaxBnd,*),nCBS(*),scffrg(*),ichfg(*),loadbf(*),
     *        loadgr(*),nstjob(*),noffg(*),indgrd(6,*),vdwrad(*),maxg(3)
c     COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      data RMC/8HMCSCF   /
c
c     Fill in some index arrays and dimension variables.
c
c     smart load balancing is based on layer one. It seems to be true that
c     it does not depend on layer anyway.
c     Next "trick" is to assume that load of a dimer is equal to the sum
c     of two monomer loads. This is not quite true if two monomers overlap
c     but should give very good idea of the dimer load. 
c
      maxl1=0
      maxnat=0
      l1fmo=0
      nefmo=0
      call viclr(maxl30,1,nfg)
      if(nbdfg.ne.0) call viclr(indatg,1,natfmo*maxabd)
      call viclr(maxg,1,3)
      ioff=1
      do ilay=1,nlayer
        do 100 ifg=1,nfg
          if(layfrg(ifg).lt.ilay) goto 100
          iifg=ifg
          iilay=ilay
          call makemol(iifg,0,0,iilay,1,0,0,0,0,0,0,.false.)
          maxl1=max(maxl1,num)
          maxnat=max(maxnat,nat)
          if(num.gt.maxl30(ifg)) maxl30(ifg)=num
          if(ilay.eq.1) numfrg(ifg)=num
          if(layfrg(ifg).eq.ilay) then
            l1fmo=l1fmo+num
            nefmo=nefmo+ne
          endif
          if(spargrid) then
            call grdbox(indgrd(1,ifg),indgrd(2,ifg),indgrd(3,ifg),
     *                  indgrd(4,ifg),indgrd(5,ifg),indgrd(6,ifg),
     *                  vdwrad,grdpad) 
            noffg(ifg)=ioff
            ioff=ioff+(indgrd(2,ifg)-indgrd(1,ifg)+1)*
     *                (indgrd(4,ifg)-indgrd(3,ifg)+1)*
     *                (indgrd(6,ifg)-indgrd(5,ifg)+1)
            maxg(1)=max(maxg(1),indgrd(2,ifg)-indgrd(1,ifg)+1)
            maxg(2)=max(maxg(2),indgrd(4,ifg)-indgrd(3,ifg)+1)
            maxg(3)=max(maxg(3),indgrd(6,ifg)-indgrd(5,ifg)+1)
          endif
  100   continue
      enddo
      nxyzg=ioff-1
      if(loadhf.eq.1) then
        if(nbody.gt.1) then
          ishift=maxl1*2
          iloop=0
          do ifg=1,nfg 
            l1i=numfrg(ifg)
            do jfg=1,ifg-1 
              iloop=iloop+1
              l1j=numfrg(jfg)
              iwrk(iloop)=l1i+l1j
              if(scffrg(ifg).eq.rmc.or.scffrg(jfg).eq.rmc)
     *          iwrk(iloop)=iwrk(iloop)+ishift
c             semidynamic jobs should be placed before all other extensions,
c             such as correlation, use the largest shift.
              if(semidyn.and.loadbf(2).ne.0.and.loadgr(2).ne.0.and.
     *           l1i+l1j.ge.loadbf(2)) then
                iwrk(iloop)=iwrk(iloop)+ishift*2
                nstjob(2)=nstjob(2)+1
              endif 
            enddo
          enddo
          call indsort(nfg2,iwrk,loadd)
c         numfrg and iwrk are destroyed after indsort
c         loadd is overwritten later in inidfmo if approximations are used.
        endif
        if(nfg3.ne.0) then
          iloop=0
          do ifg=1,nfg
            l1i=numfrg(ifg)
            do jfg=1,ifg-1
              l1ij=l1i+numfrg(jfg)
              do kfg=1,jfg-1
                iloop=iloop+1
                iwrk(iloop)=l1ij+numfrg(kfg)
              enddo
            enddo
          enddo
          call indsort(nfg3,iwrk,loadt)
        endif
        do ifg=1,nfg
          iwrk(ifg)=numfrg(ifg)
          if(scffrg(ifg).eq.rmc) iwrk(ifg)=iwrk(ifg)+maxl1
c         semidynamic jobs should be placed before all other extensions,
c         such as correlation.
          if(semidyn.and.loadbf(1).ne.0.and.loadgr(1).ne.0.and.
     *       numfrg(ifg).ge.loadbf(1)) then
            iwrk(ifg)=iwrk(ifg)+maxl1*2
            nstjob(1)=nstjob(1)+1
          endif 
        enddo
c       call indsort(nfg,numfrg,loadm)
        call indsort(nfg,iwrk,loadm)
        if(some) then
          write(6,*) 'loadm',(loadm(i),i=1,nfg)
          if(nbody.gt.1) write(6,*) 'loadd',(loadd(i),i=1,nfg2)
          if(nfg3.ne.0) write(6,*) 'loadt',(loadt(i),i=1,nfg3)
        endif
      endif
      if(some) write(6,*) 'maxl30',(maxl30(i),i=1,nfg)
c
      ichfmo=0
      docas=.false.
      do ifg=1,nfg
c       fill in max record size
        imxl30=maxl30(ifg)
        mmxl30=(imxl30*imxl30+imxl30)/2
        if(scffrg(ifg).eq.rmc) then
          docas=.true.
          mmxl30=mmxl30+imxl30*imxl30 
        else
          if(orbxch) mmxl30=imxl30*imxl30
        endif
        if(enexch) mmxl30=mmxl30+imxl30
        maxl30(ifg)=mmxl30
        ichfmo=ichfmo+ichfg(ifg)
      enddo
c
c     Subtract L1 for double counted bond-fraction joints.
c
      l1bd=0
      l0bd=0
      do ibdg=1,nbdfg
        ia=abs(iabdfg(ibdg))
        ja=abs(jabdfg(ibdg))
        ifg=indat(ia)
        jfg=indat(ja)
        ilay=min(layfrg(ifg),layfrg(jfg))
        ibdtyp=idxcao(ibdg,ilay)
        nao=nCBS(ibdtyp)
        l1bd=l1bd+nao
        l0bd=l0bd+nao-nsphel(ia,ilay)
c       write(6,*) 'wwwl00',nao,ia,ilay,nsphel(ia,ilay)
c 
c       fill in the ghost fragment array. The left end atom belongs also
c       to the right end fragment as a ghost atom.
c       Find the first empty slot.
c
        do i=maxabd,1,-1
          if(indatg(ia,i).ne.0) then
            next=i+1
            goto 200
          endif
        enddo
        next=1
  200   continue  
        if(next.gt.maxabd) then
          write(6,*) 'Increase maxabd',maxabd 
          call abrt
        endif
        indatg(ia,next)=jfg
      enddo
c
      l0fmo=l1fmo-l0bd
      l1fmo=l1fmo-l1bd
c     l0fmo is not yet finalised. spherical contaminants for each fragment
c     will be subtracted in inidfmo. Here we subtract doubly counted spherical 
c     contaminants from fractioned bonds. 
c
c     do i=1,maxabd
c       if(nbdfg.ne.0) write(6,*) (indatg(j,i),j=1,natfmo)
c     enddo 
c
c     Set charge transfer options.
c
c     do ifg=1,nfg
c       if(iand(nprfrg(ifg),8).ne.0)  ifgdon=ifg
c       if(iand(nprfrg(ifg),16).ne.0) ifgacc=ifg
c     enddo
c
      return
      end
C*MODULE fmolib  *DECK filloc
      subroutine filloc
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXSH=5000,MXGTOT=20000)
      COMMON /FMCOM / X(1)
      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 /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
c
c     fill in locfmo array and do nothing else in FMOHOP.
c     (note that this option is passed through IRSTSTP and IFMOSTP).
c     dummy arguments are passed as x and 0.
c
      call fmohop(0,0,x,x,x,x,x,x,x,x,x(liabdfg),x(ljabdfg),x(lidxCAO),
     *            x(liaglob),x,x,x,x,x,x,x,x,x(llocfmo),
     *            NSHELL,KATOM,KTYPE,KLOC,kmin,.false.)
c
      return
      end
c
C*MODULE fmolib  *DECK readcas
      subroutine readcas(nactfmo,ncasfmo)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXNORO=250)
      logical CANONC,FCORE,FORS,NOCI,EKT,LINSER
      DOUBLE PRECISION METHOD
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,NORB,
     *                NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,npreo(4)
      DATA DETWRD,DRTWRD,GENWRD/8HDET     ,8HDRT     ,8HGEN     /
      DATA GUGA,ALDET,GENCI/8HGUGA    ,8HALDET   ,8HGENCI   /
c
c     Read CAS input. Preserve the SCF value of MAXIT.
c
      MAXITsav=MAXIT
      call mcin
      IF(CISTEP.EQ.ALDET) CALL DETINP(NPFLG(1),DETWRD)
      IF(CISTEP.EQ.GUGA)  CALL DRTGEN(NPFLG(1),DRTWRD)
      IF(CISTEP.EQ.GENCI) CALL GCIINP(NPFLG(1),GENWRD)
      call fmonad(nactfmo,ncasfmo)
      MAXIT=MAXITsav
c     write(6,*) 'wwwmaxit',MAXITsav,MAXIT
c
      return
      end
c
C*MODULE fmolib  *DECK fmoord
      subroutine fmoord(v,e,iodexch,jodexch,kodexch,mapi,mapj,mapk,
     *                  enexch,l1,nai,naj,nak,iwrk1,iwrk2) 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical iodexch,jodexch,kodexch,enexch
      dimension v(l1,*),e(l1),mapi(*),mapj(*),mapk(*),iwrk1(l1),
     *          iwrk2(l1)
c
c     reorder orbitals and energies for CAS dimers/trimers.
c     iwrk2 is also used as wrk2(l1), i.e. as real array.
c     for dimers nak is zero.
c
c     call viclr(iwrk1,1,l1)
      do i=1,l1
        iwrk1(i)=0
        iwrk2(i)=i
      enddo
c     newe order 
c     first reorder I orbitals
      ind=0
      if(iodexch) ind=naj+nak
      do i=1,l1
        ii=mapi(i)
        if(ii.gt.0) then
          if(ii.le.nai) then
c           active/core
            ind=ind+1
            iwrk1(ind)=i
c         else
c           indv=indv+1
c           iwrk1(indv)=i
c         i MO is used, so throw it away from the index list.
            iwrk2(i)=0
          endif
        endif
      enddo
c     now reorder J orbitals (core or CAS)
      ind=0
      if(jodexch) ind=nai+nak
      do i=1,l1
        jj=mapj(i)
        if(jj.gt.0) then
          if(jj.le.naj) then
            ind=ind+1
            if(iwrk1(ind).ne.0) then
              write(6,*) 'Overlapping orbital indices',i,iwrk1(ind)
              call abrt
            endif
            iwrk1(ind)=i
c           i MO is used, so throw it away from the index list.
            iwrk2(i)=0
          endif
        endif
      enddo
      if(nak.ne.0) then
      ind=0
      if(kodexch) ind=nai+naj
      do i=1,l1
        kk=mapk(i)
        if(kk.gt.0) then
          if(kk.le.nak) then
            ind=ind+1
            if(iwrk1(ind).ne.0) then
              write(6,*) 'Overlapping orbital indices',i,iwrk1(ind)
              call abrt
            endif
            iwrk1(ind)=i
c           i MO is used, so throw it away from the index list.
            iwrk2(i)=0
          endif
        endif
      enddo
      endif
c     write(6,*) 'Neue Ordnung',(iwrk1(i),i=1,l1)
c     write(6,*) 'Newe Ordnung',(iwrk2(i),i=1,l1)
c     fill in virtual indices.
c     note: projected out orbitals (due to cut bonds) with lunatic energies 
c     are not systematically got rid of here. This should present no problem 
c     as one should do 
c     a dimer SCF. If one does not do that there may be a problem.
      ind=0
      do i=1,l1
        if(iwrk1(i).eq.0) then
  100     continue
          ind=ind+1
          if(iwrk2(ind).eq.0.and.ind.lt.l1) goto 100
          if(ind.eq.l1.and.i.ne.l1) call abrt
c         ran out of virtual indices?
          iwrk1(i)=ind
        endif
      enddo
c     loose ends should match
      if(ind.ne.l1) then
        write(6,*) 'Collapsed reordering',nai,naj,ind
        call abrt
      endif
c     write(6,*) 'Neue Ordnung',(iwrk1(i),i=1,l1)
      if(enexch) CALL ICOPY(L1,iwrk1,1,IWRK2,1)
      CALL REORDR(V,IWRK1,L1,L1)
      if(enexch) CALL REORDR(E,IWRK2,L1,1)
c     finally, sort separately core+active+virtuals
c     mv=nai+naj
c     nv=l1-mv
c     write(6,*) 'Sorting',mv,nv,l1
c     if(enexch) then
c       call ORDERV(v(1,   1),e(   1),iwrk1,IWRK2,34,34,l1)
c       call ORDERV(v(1,34+1),e(34+1),iwrk1,IWRK2, 2, 2,l1)
c       call ORDERV(v(1,mv+1),e(mv+1),iwrk1,IWRK2,nv,nv,l1)
c     endif
c     write(6,*) 'wwworbital-before proj'
c     call prsq(v,l1,l1,l1)
c     call projgues(l1,l2,v)
      return
      end
C*MODULE fmolib  *DECK fmonad
      subroutine fmonad(nam,nan)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical CANONC,FCORE,FORS,NOCI,EKT,LINSER
      PARAMETER (MXRT=100, MXNORO=250, MXATM=2000)
      DOUBLE PRECISION METHOD
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,NORB,
     *                NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,SDET,SZDET,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORBDT,
     *                NADET,NBDET,KDET,KSTDET,IROOT,IPURES,MAXW1,NITDET,
     *                MAXP,NCIDET,IGPDET,KSTSYM,NFTGCI
      COMMON /GUGWFN/ NFZC,NMCC,NDOC,NAOS,NBOS,NALP,NVAL,NEXT,NFZV,
     *                IFORS,IEXCIT,ICICI,NOIRR
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
c     COMMON /ORBSET/ NORBMX,NORBS,NCORBS,NLEVS,NA,NB,NC,NSYM,MSYM,
c    *                IDOCC,IVAL,IMCC,ISYM(MXAO),ICODE(MXAO),
c    *                NLCS(MXAO),LEVPT(MXAO),LEVNR(MXAO),
c    *                IOUT(MXAO),NREFS,IEXCT,NFOCI,INTACT
      data GUGA/8HGUGA    /
c
c     adjust the number of orbitals in CAS,
c     return the number of CAS active and core+active orbitals.
c
      if(cistep.eq.guga) then
        nan=NDOC+NAOS+NBOS+NALP+NVAL
        nam=NFZC+NMCC+nan
      else
        nan=NACT
        nam=NCORSV+NACT
      endif
      norb=nqmt
      return
      end
C*MODULE fmolib  *DECK fmoconv
      subroutine fmoconv(modcon0,iconfg,l1,vec,enexch,orbxch,some,tryalt
     *                  ,forcedir)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical FOCAS,SOSCF,DROPC,JACOBI,CANONC,FCORE,FORS,NOCI,EKT,LINSER
     *       ,QUD,veritas(0:1),some,enexch,orbxch,RSOSCF,rVSHIFT,tryalt,
     *        DIRSCF,FDIFF,forcedir
      dimension vec(*) 
      DOUBLE PRECISION METHOD
      PARAMETER (MXNORO=250)
      COMMON /CASOPT/ CASHFT,CASDII,NRMCAS,FOCAS,SOSCF,DROPC
      COMMON /DMPING/ SHIFTO,SHIFTV,DMPCUT,SWDIIS,DIRTHR
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /JACOBI/ JACOBI,NJAOR,ELAST,ISTAT
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,NORB,
     *                NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /QUDMC / QUDTHR,QUD
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,npreo(4)
c     COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      data veritas/.false.,.true./
c     DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF     /
c     DATA RMC/8HMCSCF   /,GVB/8HGVB     /
c
c     Set a converger for each step in FMO. Go ahead and set all SCF types.
c     Write fake initial orbitals for SOSCF (SOSCF does not need them,
c     it really needs just the density but it reads the orbitals as well).
c     tryalt: try alternative converger (that is, flip DIIS <-> SOSCF). 
c
c     SCF:
c     IF(EXTRAP) MCDEF=MCDEF+1
c     IF(DAMPH)  MCDEF=MCDEF+2
c     IF(VSHIFT) MCDEF=MCDEF+4
c     IF(RSTRCT) MCDEF=MCDEF+8
c     IF(DIIS)   MCDEF=MCDEF+16
c     IF(DEM)    MCDEF=MCDEF+32
c     IF(SOSCF)  MCDEF=MCDEF+64
c
c     MCSCF:
c     IF(FOCAS)  MCDEF=MCDEF+1024
c     IF(SOSCF)  MCDEF=MCDEF+2048
c     IF(DROPC)  MCDEF=MCDEF+4096
c     IF(CANONC) MCDEF=MCDEF+8192
c     IF(FCORE)  MCDEF=MCDEF+16384
c     IF(FORS)   MCDEF=MCDEF+32768
c     IF(NOCI)   MCDEF=MCDEF+65536
c     IF(EKT)    MCDEF=MCDEF+131072
c     IF(LINSER) MCDEF=MCDEF+262144
c     IF(JACOBI) MCDEF=MCDEF+524288
c     IF(QUD)    MCDEF=MCDEF+1048576
c
      if(modcon0.eq.-1) then
c       store the SCF defaults
        modcon0=MCONV 
        if(fdiff.and.dirscf) modcon0=logor(modcon0,256)
        if(dirscf) modcon0=logor(modcon0,512)
c       write(6,*) 'wwwaaa',modcon,dirscf
c       MCSCF options are not read in yet
      else
        modcon=modcon0
        if(iconfg.ge.0) modcon=iconfg
        MCONV=mod(modcon,1024) 
        if(tryalt) MCONV=logxor(MCONV,16+64)
c       write(6,*) 'wwwbbb',modcon,dirscf
c       write(6,*) 'new conv',MCONV
c       modcon is a global setting for all fragments.
c       iconfg is a local setting for a given fragment.
c
c       if(maswrk) write(6,*) 'Defaults',FOCAS,SOSCF,DROPC,CANONC,FCORE,
c    *                         FORS,NOCI,EKT,LINSER,JACOBI,QUD
        if(modcon.ge.1024) then
c       if less use defaults given in $MCSCF, otherwise set as told
c       note that $SCF is read in only once, whereas $MCSCF is reread for
c       each energy run, allowing not to store the default values.          
          FOCAS= veritas(mod(modcon/1024,2))
          SOSCF= veritas(mod(modcon/2048,2))
          DROPC= veritas(mod(modcon/4096,2))
          CANONC=veritas(mod(modcon/8192,2))
          FCORE= veritas(mod(modcon/16384,2))
          FORS=  veritas(mod(modcon/32768,2))
          NOCI=  veritas(mod(modcon/65536,2))
          EKT=   veritas(mod(modcon/131072,2))
          LINSER=veritas(mod(modcon/262144,2))
          JACOBI=veritas(mod(modcon/524288,2))
          QUD=   veritas(mod(modcon/1048576,2))
          if(some) write(iw,*)'Reset to',FOCAS,SOSCF,DROPC,CANONC,FCORE
     *                         ,FORS,NOCI,EKT,LINSER,JACOBI,QUD
c       else
c         SOSCF= veritas(mod(mconv/64,2))
        endif
        if(tryalt) SOSCF=.not.SOSCF
c       invert SOSCF for the alternative converger option
        rSOSCF= veritas(mod(mconv/64,2))
        rVSHIFT=veritas(mod(mconv/4,2))
        if((rSOSCF.or.rVSHIFT.or.swdiis.ne.0).and..not.orbxch) then
c         write fake initial orbitals (unit matrix) to fool SOSCF that
c         always tries to read them in. In order for this to work we must 
c         also force the first iteration to be DIIS. 
c         In some cases we should not write these dummy orbitals, e.g.
c         when doing DFT as then orbitals are essential.
          call RUNITV(l1,l1,vec)
          l3=l1*l1
          CALL dawrit(IDAF,IODA,vec,l3,15,0)
c         also write fake orbital energies, that happen to be 1,0,0,...
        endif
        if((rSOSCF.or.swdiis.ne.0).and..not.enexch) then
          call vclr(vec,1,l1)
          CALL dawrit(IDAF,IODA,vec,l1,17,0)
        endif 
        DIRSCF=iand(mconv,512).ne.0
        FDIFF=iand(mconv,256).ne.0
        if(forcedir.and.iconfg.lt.0) then
c         Do not overwrire fragment specific settings with forcedir
c         alter FDIFF only if the user has not requested DIRSCF.
          if(.not.DIRSCF) FDIFF=forcedir
          DIRSCF=forcedir
        endif
      endif
      return
      end
C*MODULE fmolib  *DECK fmorvec
      subroutine fmorvec(ifg,jfg,ijvec,naos,nmos,vec,igot)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension ijvec(2,*),vec(naos,naos)
      CHARACTER*8 cvec
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      equivalence (cvec,vecnam)
c
c     read vectors provided by the user. igot will be set to 1 if
c     orbitals were available and 0 otherwise.
c
      igot=0
      do i=1,maxvec
        if(ifg.eq.ijvec(1,i).and.jfg.eq.ijvec(2,i)) then
          WRITE(UNIT=cvec(1:8),FMT='(A5,I1)') ' $VEC',i
          call TRNRDM(IR,IW,VECNAM,naos,NMOS,VEC)
          CALL dawrit(IDAF,IODA,vec,naos*NMOS,15,0)
          if(NMOS.lt.naos) call vclr(vec(1,NMOS+1),1,naos*(naos-NMOS))
          if(maswrk) then
            if(jfg.eq.0) write(iw,*) 'Read input MOs for monomer',ifg
            if(jfg.ne.0) write(iw,*) 'Read input MOs for dimer',ifg,jfg
          endif
          igot=1
          return
        endif
      enddo
      return
      end
C*MODULE fmolib  *DECK fmosdir
      subroutine fmosdir(imode,fmodscf,fmodtrf)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical DIRSCF,DIRTRF,FDIFF,fmodscf,fmodtrf
c
c     GAMESS behaves in a bizarre way when it comes to DIRSCF and DIRTRF
c     being set differently. We want DIRSCF to be used only for RHF and
c     DIRTRF only for MCSCF and no other weird stuff.
c     imode=0: save DIRSCF and DIRTRF to fmodscf,fmodtrf; set DIRSCF,DIRTRF
c     imode=1: restore DIRSCF,DIRTRF from fmodscf,fmodtrf
c
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
c
      if(imode.eq.0) then
        fmodscf=DIRSCF
        fmodtrf=DIRTRF
c       if(scftyp.eq.RHF) DIRTRF=.false.
c       the other case (MCSCF) is harmless
      else
        DIRSCF=fmodscf
        DIRTRF=fmodtrf
      endif
      return
      end
C*MODULE fmolib  *DECK fmodor
      subroutine fmodor(enexch,VECold,vecnew,ss,ee,wrk,iwrk1,iwrk2,nact,
     *                  L1,l2,L3)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      logical enexch
      PARAMETER (ZERO=0.0D+00)
      dimension vecold(l1,l1),vecnew(l1,l1),ee(l1),ss(l2),wrk(l1),
     *          iwrk1(l1),iwrk2(l1)
c
c     find nact orbitals in vecnew having maximum overlap with VECold
c     and reorder vecnew to have the same order as VECold.
c     leave remaining orbitals unchanged
c
      call daread(IDAF,IODA,ss,L2,12,0) 
      call daread(IDAF,IODA,vecnew,L3,15,0)
c     write(6,*) 'wwwmy old'
c     call prsq(vecold,l1,l1,l1)
c     write(6,*) 'wwwmy new'
c     call prsq(vecnew,l1,l1,l1)
      call viclr(iwrk1,1,l1)
      call viclr(iwrk2,1,l1)
      ovbig=0.5D+00
      do i=1,nact 
        ii=0
        ovmax=zero
        call MTARBR(ss,l1,vecold(1,i),1,wrk,l1,1)
        do 100 j=1,l1
c         skip orbitals that were already used
          if(iwrk2(j).ne.0) goto 100
          over=abs(ddot(l1,vecnew(1,j),1,wrk,1))
          if(over.gt.ovmax) then
            ii=j
            ovmax=over 
          endif
          if(over.gt.ovbig) goto 200
  100   continue
  200   continue
c       iwrk1 keeps new orbital order. iwrk2 keeps track of used orbitals.
        iwrk1(i)=ii
        iwrk2(ii)=i
        write(6,*) 'Orbital',i,' has overlap',ovmax,' with',ii
      enddo
c     fill in remaining orbitals
      ind=nact
      do j=1,l1
        if(iwrk2(j).eq.0) then
          ind=ind+1
          iwrk1(ind)=j
        endif
      enddo 
      write(6,*) 'Neue Ordnung',(iwrk1(i),i=1,l1)
      if(enexch) call daread(IDAF,IODA,ee,L1,17,0) 
      if(enexch) CALL ICOPY(L1,iwrk1,1,IWRK2,1)
      CALL REORDR(vecnew,IWRK1,L1,L1)
      if(enexch) CALL REORDR(Ee,IWRK2,L1,1)
      if(enexch) call dawrit(IDAF,IODA,ee,L1,17,0) 
      call dawrit(IDAF,IODA,vecnew,L3,15,0) 
      return
      end
C*MODULE fmolib  *DECK matchcas
      subroutine matchcas(enexch,VECold,vecnew,ss,ee,wrk,iwrk1,iwrk2,m1,
     *                    nact,L1,l2,L3)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      logical enexch
      PARAMETER (ZERO=0.0D+00)
      dimension vecold(l1,l1),vecnew(l1,l1),ee(l1),ss(l2),wrk(l1),
     *          iwrk1(l1),iwrk2(l1)
c
c     similar to FMODOR (q.v.), but it first finds active orbitals and then
c     fills in the core.
c     M1 is the number of core+CAS active space (NACT) orbitals
c
      call daread(IDAF,IODA,ss,L2,12,0) 
      call daread(IDAF,IODA,vecnew,L3,15,0)
c     write(6,*) 'wwwmy old'
c     call prsq(vecold,l1,l1,l1)
c     write(6,*) 'wwwmy new'
c     call prsq(vecnew,l1,l1,l1)
      call viclr(iwrk1,1,l1)
      call viclr(iwrk2,1,l1)
      ovbig=0.5D+00
      do i=m1-nact+1,m1
        ii=0
        ovmax=zero
        call MTARBR(ss,l1,vecold(1,i),1,wrk,l1,1)
        do 100 j=1,l1
c         skip orbitals that were already used
          if(iwrk2(j).ne.0) goto 100
          over=abs(ddot(l1,vecnew(1,j),1,wrk,1))
          if(over.gt.ovmax) then
            ii=j
            ovmax=over 
          endif
          if(over.gt.ovbig) goto 200
  100   continue
  200   continue
c       iwrk1 keeps new orbital order. iwrk2 keeps track of used orbitals.
        iwrk1(i)=ii
        iwrk2(ii)=i
        write(6,*) 'Orbital',i,' has overlap',ovmax,' with',ii
      enddo
c     fill in core+virtual orbitals
      ind=0
      do j=1,l1
        if(iwrk2(j).eq.0) then
c         find the first empty orbital 
          do i=ind+1,l1
            ind=ind+1
            if(iwrk1(ind).eq.0) goto 500
          enddo
          write(6,*) 'CAS index could not be built',j,ind
          call abrt
  500     continue
          iwrk1(ind)=j
        endif
      enddo 
      write(6,*) 'Neue Ordnung',(iwrk1(i),i=1,l1)
      if(enexch) call daread(IDAF,IODA,ee,L1,17,0) 
      if(enexch) CALL ICOPY(L1,iwrk1,1,IWRK2,1)
      CALL REORDR(vecnew,IWRK1,L1,L1)
      if(enexch) CALL REORDR(Ee,IWRK2,L1,1)
      if(enexch) call dawrit(IDAF,IODA,ee,L1,17,0) 
      call dawrit(IDAF,IODA,vecnew,L3,15,0) 
      return
      end
C*MODULE fmolib  *DECK fmoauto
      subroutine fmoauto(nacut,indat)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension indat(*)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
c
c     automatic molecule partitioning
c
      do iat=1,natfmo
        indat(iat)=(iat-1)/nacut+1
      enddo
      return
      end
C*MODULE fmolib  *DECK fmodist
      function fmodist(ifg,jfg,kfg,lfg)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /FMCOM / X(1)
      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
c
c     Return distance between a monomer (ifg), a dimer (ifg,jfg) or
c     a trimer (ifg,jfg,kfg) and another monomer (lfg). 
c     The choice of the n-mer is accomplished by setting extra indices to 0.
c     The distance for an n-mer is defined as min(a){Ral}, a=i,j,k.
c     The distance is stored in "waals".
c     The "diagonal" distance e.g. between ifg,ifg is not defined.
c
      if(ifg.eq.lfg.or.jfg.eq.lfg.or.kfg.eq.lfg) then
c       Strictly speaking other combinations too, such as ifg and jfg,
c       but those are too unlikely by construction.
c       While mathematically it is reasonable to define R(i,i)=0 
c       it is preferred to write the program in such a way as to not
c       need it.
        write(6,*) 'Disaster in fmodist',ifg,jfg,kfg,lfg
        call abrt
      endif
      if(ifg.ge.lfg) then
        lrilfg=lrij+(ifg*ifg-3*ifg)/2+lfg
      else
        lrilfg=lrij+(lfg*lfg-3*lfg)/2+ifg
      endif
      rk=x(lrilfg)
      if(jfg.ne.0) then
        if(jfg.ge.lfg) then
          lrjlfg=lrij+(jfg*jfg-3*jfg)/2+lfg
        else 
          lrjlfg=lrij+(lfg*lfg-3*lfg)/2+jfg
        endif
        rk=min(rk,x(lrjlfg))
      endif
      if(kfg.ne.0) then
        if(kfg.ge.lfg) then
          lrklfg=lrij+(kfg*kfg-3*kfg)/2+lfg
        else
          lrklfg=lrij+(lfg*lfg-3*lfg)/2+kfg
        endif
        rk=min(rk,x(lrklfg))
      endif
      fmodist=rk
      return
      end
C*MODULE fmolib  *DECK fmodist3
      subroutine fmodist3(ifg,jfg,kfg,rmin,rmax)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /FMCOM / X(1)
      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
c
c     Return inter-fragment distances for a trimer:
c     rd the distance between two closest monomers, forming a dimer,
c     rm the distance between this dimer and the remaining monomer (indm).
c     The distance is stored in "waals".
c
      if(ifg.ge.jfg) then
        lrijfg=lrij+(ifg*ifg-3*ifg)/2+jfg
      else
        lrijfg=lrij+(jfg*jfg-3*jfg)/2+ifg
      endif
      rij=x(lrijfg)
      if(ifg.ge.kfg) then
        lrikfg=lrij+(ifg*ifg-3*ifg)/2+kfg
      else
        lrikfg=lrij+(kfg*kfg-3*kfg)/2+ifg
      endif
      rik=x(lrikfg)
      if(jfg.ge.kfg) then
        lrjkfg=lrij+(jfg*jfg-3*jfg)/2+kfg
      else
        lrjkfg=lrij+(kfg*kfg-3*kfg)/2+jfg
      endif
      rjk=x(lrjkfg)
        rmin=rjk
        rmax=min(rij,rik)
      if(rij.lt.rik.and.rij.lt.rjk) then
        rmin=rij
        rmax=min(rik,rjk)
      else if(rik.lt.rij.and.rik.lt.rjk) then
        rmin=rik
        rmax=min(rij,rjk)
      endif
      return
      end
C*MODULE fmolib  *DECK fmosdist
      function fmosdist(iifg0,jjfg0,iifgg,jjfgg,ifg,jfg,kfg,lfg)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /FMCOM / X(1)
      parameter (huge=1.0D+32)
      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 /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      dimension iifgg(*),jjfgg(*)
c
c     Return distance between a monomer (ifg), a dimer (ifg,jfg) or
c     a trimer (ifg,jfg,kfg) and another monomer (lfg).
c     The choice of the n-mer is accomplished by setting extra indices to 0.
c     This subroutine defines distance as a function of iifg,jjfg that are
c     row and column fragments of the ESP Vi,j matrix, where iifg and jjfg are
c     fragments to which i and j belong. In other words, in case of a dimer
c     there are 3 blocks in the triangular matrix Vij:
c     ii block where i,j are from ifg. The distance fmosdist is then ril.
c     ij block where i,j are from ifg,jfg. The distance is then min(ril,rjl).
c     jj block where i,j are from jfg. The distance fmosdist is then rjl.
c     Note that some i (or j) may belong to more than one fragment (on the
c     cloven bonds; e.g. to both ifg and jfg). That does not pose a problem 
c     because in this case rij (the distance between ifg and jfg) is zero 
c     and ril=rjl, so one can use either ril or rjl (that is, assume that i
c     belongs to only one fragment).
c     Trouble occurs for ghost atoms: they belong (according to indat) to 
c     neither of (ifg,jfg,kfg). The solution is to use ghost index arrays
c     that contain alternative fragment IDs. 
c     fmosdist produces the same answer as fmodist only for monomers (jfg=kfg=0)
c     and in all cases fmosdist>=fmodist with the same arguments ifg,jfg,kfg,lfg
c     The distance is stored in "waals".
c
c     First redirect ghost atoms to one of ifg,jfg,kfg.
c
      fmosdist=huge 
      maxabd0=maxabd
      if(nbdfg.eq.0) maxabd0=0
      do 1000 i=0,maxabd0
        iifg=iifg0
        if(i.gt.0) then
          iifgi=iifgg(i)
          if(iifgi.eq.0) goto 100
          if(iifgi.eq.ifg) then
            iifg=ifg
            goto 200
          endif
          if(iifgi.eq.jfg) then
            iifg=jfg
            goto 200
          endif
          if(iifgi.eq.kfg) then
            iifg=kfg
            goto 200
          endif
c         if not found loop again 
          goto 1000
  200     continue
        endif
        do 2000 j=0,maxabd0
          jjfg=jjfg0
          if(j.gt.0) then
            jjfgi=jjfgg(j)
            if(jjfgi.eq.0) goto 1000
            if(jjfgi.eq.ifg) then
              jjfg=ifg
              goto 300
            endif
            if(jjfgi.eq.jfg) then
              jjfg=jfg
              goto 300
            endif
            if(jjfgi.eq.kfg) then
              jjfg=kfg
              goto 300
            endif
c           if not found loop again
            goto 2000
  300       continue
          endif
c
c     Now do the real distance work. 
c
      if(iifg.eq.ifg.or.jjfg.eq.ifg) then 
        if(ifg.ge.lfg) then
          lrilfg=lrij+(ifg*ifg-3*ifg)/2+lfg
        else
          lrilfg=lrij+(lfg*lfg-3*lfg)/2+ifg
        endif
        ril=x(lrilfg)
        if(fmosdist.gt.ril) fmosdist=ril
      endif
c     with sane arguments if jfg is 0 (monomer) iifg and jjfg are never zeros
      if(iifg.eq.jfg.or.jjfg.eq.jfg) then
        if(jfg.ge.lfg) then
          lrjlfg=lrij+(jfg*jfg-3*jfg)/2+lfg
        else
          lrjlfg=lrij+(lfg*lfg-3*lfg)/2+jfg
        endif
        rjl=x(lrjlfg)
        if(fmosdist.gt.rjl) fmosdist=rjl
      endif
c     with sane arguments if kfg is 0 (1,2-mer) iifg and jjfg are never zeros
      if(iifg.eq.kfg.or.jjfg.eq.kfg) then
        if(kfg.ge.lfg) then
          lrklfg=lrij+(kfg*kfg-3*kfg)/2+lfg
        else
          lrklfg=lrij+(lfg*lfg-3*lfg)/2+kfg
        endif
        rkl=x(lrklfg)
        if(fmosdist.gt.rkl) fmosdist=rkl
      endif
 2000 continue 
 1000 continue 
  100 continue
      if(fmosdist.eq.huge) call abrt
c     abort means indices (calling arguments) are out of match:
c     iifg and jjfg must be equal to some of (ifg,jfg,kfg).
      return
      end
C*MODULE fmolib  *DECK monbsr
      subroutine monbsr(nat,ich,mul,num,nqmt,ne,na,nb,nshell,ngau,enucr,
     *                  nat0,ich0,mul0,num0,nqmt0,ne0,na0,nb0,
     *                  nshell0,ngau0,enucr0)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c
c     monomer basis save / restore
c
      nat0=nat
      ich0=ich
      mul0=mul
      num0=num
      nqmt0=nqmt
      ne0=ne
      na0=na
      nb0=nb
      nshell0=nshell
      ngau0=ngau
      enucr0=enucr
c
      return
      end
C*MODULE fmolib  *DECK fmoprr
      subroutine fmoprr(indat,iabdfg,jabdfg,fmozan,fmoc,vdwr,rij,prall,
     *                  prtdst,nerr)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension indat(*),iabdfg(*),jabdfg(*),fmozan(*),fmoc(3,*),vdwr(*)
     *         ,rij(*),prtdst(3)
      logical prall
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
C
C     ----- print interfragment DISTANCES -----
C
      nerr=0
      if(nfg.eq.1) return
      if(prall) then 
        WRITE (IW,9000)
        INCR = 8
        IF(NPRINT.EQ.6) INCR=16
        do i=nfg,2,-incr
          ni=min(incr,i-1)
          write(iw,9005) (i-k,k=0,ni-1)
          do j=1,i-1
            nk=min(incr,i-j)-1
            write(iw,9010) j,(rij(((i-k)*(i-k)-3*(i-k))/2+j+1),k=0,nk)
          enddo
        enddo
        write(6,9020)
      endif
      if(prtdst(1).ne.0.or.prtdst(2).ne.0.or.prtdst(3).ge.0) then
        write(6,9100)
        loop=0 
        do i=2,nfg
          do j=1,i-1
            loop=loop+1
            rr=rij(loop)
            if(rr.lt.prtdst(2).and.rr.ne.0) then
              write(iw,9110) i,j,rr
            else if(rr.lt.prtdst(1)) then
              write(iw,9120) i,j,rr
            endif
c
c           Count fractioned bonds for a pair of fragments.
c
            ibdfg1=0
            nfract=0
            do ibdfg=1,nbdfg
              ifg=indat(abs(iabdfg(ibdfg)))
              jfg=indat(abs(jabdfg(ibdfg)))
              if(ifg.eq.i.and.jfg.eq.j.or.ifg.eq.j.and.jfg.eq.i) then
c               output cross-links if more than one is present.
                if(nfract.eq.1) write(iw,9210) i,j,nfract,ibdfg1,
     *                           abs(iabdfg(ibdfg1)),abs(jabdfg(ibdfg1))
                nfract=nfract+1
                ibdfg1=ibdfg
                if(nfract.ge.2) write(iw,9210) i,j,nfract,ibdfg,
     *                           abs(iabdfg(ibdfg)),abs(jabdfg(ibdfg))
              endif
            enddo
c
c           Check fragmentation carefully. 
c
            if(rr.le.prtdst(3).and.nfract.eq.0) then
c             In the majority of cases nfract is 1.
c             nfract>1 is a weird case that may be permissible sometimes.
c             E.g. two sulphur-bridged fragments or otherwise twice 
c             connected pair of fragments. 
c             nfract=0 can be a major mistake in $FMOBND.
c             There are hydrogen bonds etc that may justify nfract=0.
c             This is a suspicious fragment pair:close but no bond in $FMOBND.
c             Consider only heavy (non-H) atoms (no ghosts).
c             It may be a good idea to check hydrogens too (for steric
c             hindrances due to poor modelling) but here we only check bonds
c             that have to be but are not fractioned. 
              rijh=1.0D+30
              imin=0
              jmin=0
              do iat=1,natfmo
                ifg=indat(iat) 
                if(ifg.eq.i) then
                  ian=int(fmozan(iat)+1.0D-03)
                  if(ian.gt.1) then
                    ri=vdwr(ian)
                    xi=fmoc(1,iat)
                    yi=fmoc(2,iat)
                    zi=fmoc(3,iat)
                    do jat=1,natfmo
                      jfg=indat(jat)
                      if(jfg.eq.j) then
                        jan=int(fmozan(jat)+1.0D-03)
                        if(jan.gt.1) then
                          rj=vdwr(jan)
                          r12=sqrt((xi-fmoc(1,jat))**2+
     *                             (yi-fmoc(2,jat))**2+
     *                             (zi-fmoc(3,jat))**2)/(ri+rj)
                          if(r12.lt.rijh) then
                            rijh=r12
                            imin=iat 
                            jmin=jat 
                          endif
                        endif
                      endif
                    enddo
                  endif
                endif
              enddo
c             rijh is now the shortest distance between heavy atoms in I,J.
              if(rijh.le.prtdst(3)) then
                write(iw,9200) i,j,rr,rijh,imin,jmin
                nerr=nerr+1
              endif
c           else if(nfract.gt.1) then
c             write(iw,9210) i,j,rr,nfract
            endif
          enddo
        enddo
c       if(prtdst(3).ne.0) call timit(1)
      endif
      RETURN
 9000 FORMAT(/,10X,'INTERFRAGMENT DISTANCES',/,10X,30(1H-))
 9005 FORMAT(/1x,5x,10(I5,3x))
 9010 FORMAT(1x,I5,10F8.3)
 9020 FORMAT(1x)
 9100 FORMAT(/1x,'Connected fragment pairs, distance relative to vdW',
     *           ' radii',/)
 9110 FORMAT(1x,2I6,F10.5,' very close!')
 9120 FORMAT(1x,2I6,F10.5)
 9200 FORMAT(1x,'Fragments',2I6,' are separated by',F7.3,' (all atoms)',
     *          F7.3,' (heavy only).',/1x,'Atoms',2I6,' may have to ',
     *          'have a bond between them defined in $FMOBND.',/)
 9210 FORMAT(1x,'Warning: frgs',2I6,' are cross-linked(',I1,'): bond',
     *          I6,' atoms',2I7)
      end
C*MODULE fmolib  *DECK setindbd
      subroutine setindbd(iabdfg,jabdfg,indbd)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension iabdfg(*),jabdfg(*),indbd(maxabd,*)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
c
c       fill in broken bond index for each atom (there can be several broken
c       bonds for a given atom (if none, store 0).
c
      call viclr(indbd,1,natfmo*maxabd)
      do ibdg=1,nbdfg
        ia=abs(iabdfg(ibdg))
        ja=abs(jabdfg(ibdg))
        ii=1
        jj=1
        do k=1,maxabd
          if(indbd(k,ia).ne.0) ii=k+1
          if(indbd(k,ja).ne.0) jj=k+1
        enddo
        if(ii.gt.maxabd.or.jj.gt.maxabd) then
          write(6,*) 'too many abds',ii,jj,maxabd
          call abrt
        endif
        indbd(ii,ia)=ibdg 
        indbd(jj,ja)=ibdg
      enddo
c     write(6,6666) ((indbd(k,ia),k=1,maxabd),ia=1,natfmo)
c6666 format(1x,'indbd ',4I6)
      return
      end
C*MODULE fmolib  *DECK fmogind
      subroutine fmogind(maxi,indat,indatg)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension indat(*),indatg(*)
      logical GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
c
c     process Gaussian-like INDAT, indicated by indat(1)=0 (therefore, skip
c     indat(1))
c
      ifg=1
      nifg=0
      natot=0
      i=1
  100 continue
        i=i+1
        if(i.ge.maxi) goto 200
        now=indat(i)
        if(now.eq.0) then
          if(ifg.eq.nfg) goto 200
          if(nifg.eq.0) then
            if(maswrk) write(iw,*) 'No atoms in fragment',ifg
            call abrt
          endif
          ifg=ifg+1
          natot=natot+nifg
          nifg=0
          goto 100
        endif
        if(indat(i+1).lt.0) then
          i=i+1
          next=abs(indat(i))
        else
          next=now
        endif
        do j=now,next
          nifg=nifg+1
          indatg(j)=ifg
        enddo
      goto 100
  200 continue
      natot=natot+nifg
      if(ifg.ne.nfg.or.natot.ne.natfmo) then
        if(maswrk) write(iw,9000) ifg,nfg,natot,natfmo
        call abrt
      endif
      call icopy(natfmo,indatg,1,indat,1)
      return
 9000 format(/1x,'Bad indat: nfg(read,expected)=',2I5,' natfmo=',2I5,/)
      end
C*MODULE fmolib  *DECK makefg
      subroutine makefg(indat,iabdfg,jabdfg,fmozan,fmoc,ialoc,iatfrg,
     *                 indfrg,indgfrg,natfrg,nat0frg,ianfrg,zanfrg,cfrg)
c
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXSH=5000,MXATM=2000,one=1.0D+00)
      dimension indat(*),iabdfg(*),jabdfg(*),fmozan(*),fmoc(3,*),
     *          ialoc(*),iatfrg(*),indfrg(*),indgfrg(*),natfrg(*),
     *          nat0frg(*),ianfrg(*),zanfrg(*),cfrg(3,*)
      logical GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      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 /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
c
c     Prepare fragment specific data. This subroutine is not parallelised!
c
      if(maswrk) write(iw,*) 'Initialising fragments' 
      call timit(1)
c
c     first check index.
c
      do i=1,natfmo
        if(indat(i).gt.nfg.or.indat(i).le.0) then
          if(maswrk) write(6,*) 'Wrong indat:',i,indat(i),nfg
          call abrt
        endif
      enddo
c
      natfmob=natfmo+nbdfg
      ind=0
      indp=0
      do ifg=1,nfg
        indp=ind
        indfrg(ifg)=ind+1
        ifg16=ishft(ifg,16)
c       Find normal atoms that belong to fragment IFG 
        do i=1,natfmo
          if(indat(i).eq.ifg) then
            ind=ind+1
            iatfrg(ind)=i
            cfrg(1,ind)=fmoc(1,i) 
            cfrg(2,ind)=fmoc(2,i) 
            cfrg(3,ind)=fmoc(3,i) 
            zanfrg(ind)=fmozan(i)
            ianfrg(ind)=int(fmozan(i)+0.5D+00)
            ialoc(i)=ind-indp+ifg16
          endif
        enddo
        indgfrg(ifg)=ind+1
        nat0frg(ifg)=ind-indp
        if(nat0frg(ifg).eq.0) then
          write(6,*) 'Fragment',ifg,' has no atoms assigned!'
          call abrt
        endif
        do ibdfg=1,nbdfg
          ia=abs(iabdfg(ibdfg))
          ja=jabdfg(ibdfg)
c         Add ghost atoms that belong to fragment IFG 
          if(indat(ja).eq.ifg) then
            ind=ind+1
            iatfrg(ind)=ia
            cfrg(1,ind)=fmoc(1,ia)
            cfrg(2,ind)=fmoc(2,ia)
            cfrg(3,ind)=fmoc(3,ia)
            zanfrg(ind)=one
            ianfrg(ind)=int(fmozan(ia)+0.5D+00)
            ialoc(natfmo+ibdfg)=ind-indp+ifg16
          endif
c         Subtract the ghost atom charges from fragment IFG 
          if(indat(ia).eq.ifg) then
            ind1=indp+iand(ialoc(ia),65535)
            zanfrg(ind1)=zanfrg(ind1)-one
          endif
        enddo
        natfrg(ifg)=ind-indp
      enddo
      if(ind.ne.natfmob) then
        write(iw,*) 'Invalid INDAT:',ind,natfmob
        call abrt
      endif
c     create atom mapping for C1 
      do i=1,MXATM
        mapctr(i,1)=i
      enddo 
      do i=1,MXSH
        MAPSHL(i,1)=i
      enddo 
      call timit(1)
      return
      end
C*MODULE fmolib  *DECK addfrg
      subroutine addfrg(nati,indi,zanfrg,cfrg,ianfrg,iatfrg,iaglob)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      dimension iaglob(*),zanfrg(*),cfrg(3,*),ianfrg(*),iatfrg(*)
c
      if(nat+nati.gt.MXATM) then
        if(maswrk) write(iw,9100) nat+nati,MXATM
        call abrt
      endif
      ii=nat+1 
      call dcopy(3*nati,cfrg(1,indi),1,c(1,ii),1)
      call dcopy(nati,zanfrg(indi),1,zan(ii),1)
      call icopy(nati,ianfrg(indi),1,ian(ii),1)
      call icopy(nati,iatfrg(indi),1,iaglob(ii),1)
      nat=nat+nati 
      return
 9100 format(/1x,'GAMESS must be recompiled with a larger MXATM',2I7/)
      end
C*MODULE fmolib  *DECK addgho
      subroutine addgho(natg,zang,cg,iang,indg,natj,indj,natjp,natk,indk
     *                 ,natkp,iaglob)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      dimension zang(*),cg(3,*),iang(*),indg(*),indj(*),indk(*),
     *          iaglob(*)
c     
      do i=1,natg
        iat=indg(i)
        zi=zang(i)
        ifound=0
        do j=1,natj
          if(iat.eq.indj(j)) then
            ifound=j+natjp
            goto 100
          endif
        enddo
        do k=1,natk
          if(iat.eq.indk(k)) then
            ifound=k+natkp
            goto 100
          endif
        enddo
  100   continue
        if(ifound.ne.0) then
          zan(ifound)=zan(ifound)+zi
        else
          nat=nat+1
          if(nat.gt.MXATM) then
            if(maswrk) write(iw,9100) nat,MXATM
            call abrt
          endif
          zan(nat)=zi
          c(1,nat)=cg(1,i)
          c(2,nat)=cg(2,i)
          c(3,nat)=cg(3,i)
          ian(nat)=iang(i)
          iaglob(nat)=iat
        endif
      enddo
      return
 9100 format(/1x,'GAMESS must be recompiled with a larger MXATM',2I7/)
      end
C*MODULE fmolib  *DECK moldim
      subroutine moldim(molfrg,iwrk)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension molfrg(nfg),iwrk(nfg)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
c
c     expand molfrg from compact to indexed format:
c     that is replace molfrg(1)=3,5,0 by
c     molfrg(1)=0,0,1,0,1,(0,).
c
      ires=0
      if(molfrg(1).eq.0) ires=-1 
      do ifg=1,nfg
        iwrk(ifg)=ires
      enddo
      if(ires.eq.0) then 
        do ifg=1,nfg
          mfg=molfrg(ifg)
          if(mfg.ne.0) iwrk(mfg)=1
        enddo
      endif
      call icopy(nfg,iwrk,1,molfrg,1)
      return
      end
C*MODULE fmolib  *DECK adjconv
      subroutine adjconv(iter,denmax,itol0,icut0,CONVHF0,fmajor)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical GOPARR,DSKWRK,MASWRK,fmajor
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,npreo(4)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      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
c
c     adjust convergence SCF parameters depending on monomer SCF convergence
c
      cnvmin=1.0d-04
      icutmin=6
      if(denmax.eq.0) then
c       save original convergence values
        icut0=icut
        itol0=itol
        CONVHF0=CONVHF
c       set fairly tight values
        CONVHF=1.0d-05
        icut=8
        itol=16
      else 
        if(fmajor) then
c         Force full values.
          icut=icut0
          itol=itol0
          CONVHF=CONVHF0
        else
c       determine new SCF convergence
c       Do not permit it to be ridiculous
        CONVH=min(cnvmin,denmax/cnvdmp)
c       nor go beyond the desired value
        CONVH=max(CONVH,CONVHF0) 
c       nor increase relative to previous value: otherwise weird oscillations
c       may occur.
        CONVHF=min(CONVHF,CONVH)
c       icutn=max(icutmin,int(-log10(CONVHF))+2)
        icutn=max(icutmin,icut0-int(log10(denmax/convfg)))
c       2 is some safeguard, set integral accuracy to be 100 times better 
c       than SCF
        icutn=min(icutn,icut0)
        icut=max(icut,icutn)
c       itol is more difficult to set, try twice icut.
        itol=min(icut*2,itol0)
        endif
      endif
      if(maswrk) write(iw,9000) iter,CONVHF,icut,itol
      return
 9000 format(1x,'Resetting SCF convergence for iter',I3,' to CONV=',
     *          E8.2,', ICUT=',I2,', ITOL=',I3)
      END
C*MODULE fmolib  *DECK fmolag
      subroutine fmolag(eps,V,wrk,L1,l2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (half=0.5D+00)  
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /FMCOM / X(1)
      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
      dimension eps(l2),v(l1,l1),wrk(l1)
c
c     Correct the Lagrangian subtracting diagonal ESP elements in MO basis 
c     (it looks like adding because eps has a minus in front).
c     lfmoespb is used as temporary storage L2.
c     lfmoespa contains ESP (as stored in FMOESP).
c     Lagrangian correctionis 1/2*DVD.
c     At present the hermiticity of D is not used.
c
c       write(6,*) '1e+2e ESP in grd is',l1,na
c       call prtril(x(lfmoespa),l1)
c       write(6,*) 'orbs'
c       call prsq(v,l1,l1,l1)
c     call TFDIAG(x(lfmoespa),V,wrk,WRK1,L1,L2,l1,na)
      CALL daread(IDAF,IODA,x(lfmoespb),l2,16,0)
      call CPYTSQ(x(lfmoespb),V,L1,1)
      call TFTRI(x(lfmoespb),x(lfmoespa),v,WRK,l1,l1,l1)
      call daxpy(l2,half,x(lfmoespb),1,eps,1)
c       call prtri(wrk,na)
      return
      END
C*MODULE fmolib  *DECK fndcntr
      SUBROUTINE fndcntr(fmozan,fmoc,iwhere,Xc,yc,zc)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter (zero=0.0D+00)
      dimension fmozan(*),fmoc(3,*),cc(3)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
c
c     Find the centre of mass/charge as needed for ES moments.
c     The centre of iwhere=0 mass, 1 charge, -1 none.
c     At present no parallelisation. 
c
      call vclr(cc,1,3)
      centre=zero
      if(iwhere.eq.1) then
        do iatm=1,natfmo 
          do i=1,3
            cc(i)=cC(i)+fmoC(i,iatm)*fmoZAN(iatm)
          enddo
          centre=centre+fmoZAN(iatm)
        enddo
      endif
      if(iwhere.ge.0) then
        if(abs(centre).lt.1.0D-08) then
          xc=zero
          yc=zero
          zc=zero
        else
          xc=cc(1)/centre
          yc=cc(2)/centre
          zc=cc(3)/centre
        endif
      endif
      RETURN
      END
C*MODULE fmolib  *DECK fmossdist
      function fmossdist(iifg0,jjfg0,iifgg,jjfgg,ifg,jfg,kfg,lfg,bimer)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (huge=1.0D+32)
      logical bimer(3),needi,needj,needk 
      COMMON /FMCOM / X(1)
      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 /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      dimension iifgg(*),jjfgg(*)
c
c     Return distance between a monomer (ifg), a dimer (ifg,jfg) or
c     a trimer (ifg,jfg,kfg) and another monomer (lfg).
c     The choice of the n-mer is accomplished by setting extra indices to 0.
c     This subroutine defines distance as a function of iifg,jjfg that are
c     row and column fragments of the ESP Vi,j matrix, where iifg and jjfg are
c     fragments to which i and j belong. In other words, in case of a dimer
c     there are 3 blocks in the triangular matrix Vij:
c     ii block where i,j are from ifg. The distance fmosdist is then ril.
c     ij block where i,j are from ifg,jfg. The distance is then min(ril,rjl).
c     jj block where i,j are from jfg. The distance fmosdist is then rjl.
c     Note that some i (or j) may belong to more than one fragment (on the
c     cloven bonds; e.g. to both ifg and jfg). That does not pose a problem 
c     because in this case rij (the distance between ifg and jfg) is zero 
c     and ril=rjl, so one can use either ril or rjl (that is, assume that i
c     belongs to only one fragment).
c     Trouble occurs for ghost atoms: they belong (according to indat) to 
c     neither of (ifg,jfg,kfg). The solution is to use ghost index arrays
c     that contain alternative fragment IDs. 
c     fmosdist produces the same answer as fmodist only for monomers (jfg=kfg=0)
c     and in all cases fmosdist>=fmodist with the same arguments ifg,jfg,kfg,lfg
c     The distance is stored in "waals".
c
c     Note on bimer: for connected trimers, all bimers must be set to true 
c     (not just the two directly connected).
c
c     First redirect ghost atoms to one of ifg,jfg,kfg.
c
      iifg=iifg0
      jjfg=jjfg0
      if(nbdfg.ne.0) then
        if(iifg0.ne.ifg.and.iifg0.ne.jfg.and.iifg0.ne.kfg) then
c         iifg0 is a ghost atom to be redirected
          do i=1,maxabd
            iifgi=iifgg(i)
            if(iifgi.eq.0) goto 100
            if(iifgi.eq.ifg) then
              iifg=ifg
              goto 100
            endif
            if(iifgi.eq.jfg) then
              iifg=jfg
              goto 100
            endif
            if(iifgi.eq.kfg) then
              iifg=kfg
              goto 100
            endif
          enddo
  100     continue
        endif
        if(jjfg0.ne.ifg.and.jjfg0.ne.jfg.and.jjfg0.ne.kfg) then
c         jjfg0 is a ghost atom to be redirected
          do i=1,maxabd
            jjfgi=jjfgg(i)
            if(jjfgi.eq.0) goto 200
            if(jjfgi.eq.ifg) then
              jjfg=ifg
              goto 200
            endif
            if(jjfgi.eq.jfg) then
              jjfg=jfg
              goto 200
            endif
            if(jjfgi.eq.kfg) then
              jjfg=kfg
              goto 200
            endif
          enddo
  200     continue
        endif
c
c     Now force "bound fragments". This means if rij=0 then i and j become one,
c     so if i fragment is used to get distance, j is used two, irrespective o
c     what iifg and jjfg directly imply (they both can point just to i). 
c
      endif
c
c     Now do the real distance work. 
c
      fmossdist=huge 
      needi=iifg.eq.ifg.or.jjfg.eq.ifg
      needj=iifg.eq.jfg.or.jjfg.eq.jfg
      needk=iifg.eq.kfg.or.jjfg.eq.kfg
      if(needi.or.needj.and.bimer(1).or.needk.and.bimer(2)) then
        if(ifg.ge.lfg) then
          lrilfg=lrij+(ifg*ifg-3*ifg)/2+lfg
        else
          lrilfg=lrij+(lfg*lfg-3*lfg)/2+ifg
        endif
        fmossdist=x(lrilfg)
      endif
c     with sane arguments if jfg is 0 (monomer) iifg and jjfg are never zeros
      if(needj.or.needi.and.bimer(1).or.needk.and.bimer(3)) then
        if(jfg.ge.lfg) then
          lrjlfg=lrij+(jfg*jfg-3*jfg)/2+lfg
        else
          lrjlfg=lrij+(lfg*lfg-3*lfg)/2+jfg
        endif
        rjl=x(lrjlfg)
        if(fmossdist.gt.rjl) fmossdist=rjl
      endif
c     with sane arguments if kfg is 0 (1,2-mer) iifg and jjfg are never zeros
      if(needk.or.needi.and.bimer(2).or.needj.and.bimer(3)) then
        if(kfg.ge.lfg) then
          lrklfg=lrij+(kfg*kfg-3*kfg)/2+lfg
        else
          lrklfg=lrij+(lfg*lfg-3*lfg)/2+kfg
        endif
        rkl=x(lrklfg)
        if(fmossdist.gt.rkl) fmossdist=rkl
      endif
      if(fmossdist.eq.huge) call abrt
c     abort means indices (calling arguments) are out of match:
c     iifg and jjfg must be equal to some of (ifg,jfg,kfg).
      return
      end
C*MODULE fmolib  *DECK fndchr
      function ifndchr(string,len,chr)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      character*1 string(len),chr
c
c     find chr in string and return its position (0 if not found)
c
      do i=1,len
        if(string(i).eq.chr) then
          ifndchr=i
          return
        endif
      enddo
      ifndchr=0
      RETURN
      END
C*MODULE fmolib  *DECK convlmo
      subroutine convlmo(nao,nlmo,nlmo0,in0,iaprjo,japrjo,CoreAO,libish)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (shiftb=1.0D+06,tolb=1.0D+01)
      PARAMETER (MXSH=5000,MXGTOT=20000)
      PARAMETER (MXATM=2000,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 /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 /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /FMCOM / X(1)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      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 iaprjo(MaxCAO),japrjo(MaxCAO),CoreAO(MaxCBS,MaxCAO,*),
     *          nlmo0(0:1),libish(MAXNZ,maxbas,*)
c
c     convert LMOs.   
c     for now assume 1 layer, 1 basis set
c     shiftb in this subroutine may be different from FMOHOP and such.
c     It is only used to tell black from white, so any large value will do. 
c
      if(maxbas.ne.1.or.nlayer.ne.1) call abrt
c
      nao2=(nao*nao*nao)/2
      l2=(NUM*NUM+NUM)/2
      CALL VALFM(LOADFM)
      lss=LOADFM+1
      ls=lss+l2
      lsq=ls+nao2
      ld=lsq+nao*nao
      lp=ld+nao*nao
      lq=lp+nao2
      le=lq+nao*nao
      lscr=le+nao
      liwrk=lscr+nao*8
      last=liwrk+nao
      NEED = LAST- LOADFM -1
      CALL GETFM(NEED)
      CALL daread(IDAF,IODA,x(lss),l2,12,0)
c
c     Find what atom that is. For now hardwire to C.
c 
      iz=6
      ilay=1
      ibas=1
      jz=0
      do jat=1,nat
        jz=int(fzan(jat)+0.5D+00)
        jlay=mod(llay(jat),100)
        jbas=llay(jat)/100
        if(jz.eq.iz.and.ilay.eq.jlay.and.ibas.eq.jbas) goto 100
      enddo  
  100 continue
      if(jz.eq.0) then
        write(6,*) 'Bonding atom not found'
        call abrt
      endif
      ii=libish(jz,ibas,ilay)
cnb   for now use kloc since lloc is not saved
      iloc=kloc(ii)
c     copy a triangular block from the triangular matrix SS
      iloop=(iloc*iloc+iloc)/2-1
      jloop=0
      do i=1,nao
        call dcopy(i,x(lss+iloop),1,x(ls+jloop),1)
        iloop=iloop+iloc+i-1
        jloop=jloop+i
      enddo
c     write(6,*) 'wwwaha',NAT,nao,nlmo,ii,iloc
c     call prtri(x(lss),num)
c     call prtri(x(ls),nao)
c
c     Compute the Q matrix for one atom.
c     No d functions yet!
c
      call dcopy(nao2,x(ls),1,x(lp),1) 
      nao0=nao
      call QMATRX(x(lp),x(lq),x(le),x(lscr),x(lIWrk),nbo0,nao0,nao0,
     *            .false.)
      call MTARBR(x(ls),nao,x(lq),nao0,x(lsq),nao,1)
c
c     Loop over sides: 0 left, 1 right
c
      do iside=0,1
c
c       Form B=C*b*Ct
c       where C are LMO LCAO coefficients and b is the diagonal matrix of B 
c       constants.
c
        call vclr(x(ld),1,nao2)
        nmos=0
        do imo=1,nlmo
          if(iside.eq.0.and.iaprjo(imo).ne.0 .or.
     *       iside.ne.0.and.japrjo(imo).ne.0) then
            iloop=0
            do i=1,nao
              call daxpy(i,shiftb*CoreAO(i,imo,in0),CoreAO(1,imo,in0),1,
     *                   x(ld+iloop),1)
              iloop=iloop+i
            enddo
            nmos=nmos+1
          endif
        enddo
c       call prtri(x(ld),nao)
c       Form Pq=Qt * P * Q = Qt * St * B * S * Q
        call TFTRI(x(lp),x(ld),x(lsq),x(liwrk),nao0,nao,nao)
c       diagonalise Pq
        nkept=nao0-nmos
        IGERR = 0
        CALL GLDIAG(nao,nao0,nao0,x(lp),x(lSCR),x(lE),x(ld),
     *              IGERR,x(lIWRK))
        IF (IGERR .NE. 0) CALL ABRT
        CALL TFSQB(X(Ld),X(LQ),X(LSCR),nao0,nao,nao)
c       check upper eigenvalues
        do imo=nao0-nmos+1,nao0
          if(x(lE+imo-1).lt.shiftb/tolb) then
            write(6,*) 'P eigenvalue too small?',x(lE+imo-1),imo
            call abrt 
          endif
        enddo
c       write(6,*) 'wwwZ',iside
c       call prsq(X(Ld),nao0,nao,nao)
        ishift=nlmo0(iside)
        do i=1,nao0
          call dcopy(nao,x(ld+(i-1)*nao),1,CoreAO(1,i,ishift),1)
        enddo
c       call dcopy(MaxCBS*nao,x(ld),1,CoreAO(1,1,ishift),1)
        nlmo0(iside)=nkept
        write(6,*) 'Generated ',nkept,'orbitals for side',iside
      enddo
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE fmolib  *DECK dolend
      logical function dolend(str)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      character*99 str
      character*5 str0
c
c     check if str is equal to $end (ignoring case and following characters).
c 
      do i=1,5
        str0(i:i)=str(i:i) 
      enddo
      call UPRCAS(str0,5)
      dolend=str0.eq.' $END' 
      RETURN
      END
c
C*MODULE fmolib  *DECK glolat
      SUBROUTINE glolat(loctat,libish,libnsh,izbas,fmozan)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MAXL=5,MAXNZ=137,MaxLay=5,
     *           MXSFMO=MaxLay*40,MXGFMO=MaxLay*100,MXAFMO=MaxLay*10)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      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 loctat(*),libish(MAXNZ,maxbas,*),libnsh(MAXNZ,maxbas,*),
     *          izbas(*),fmozan(*) 
c
c     Fill the array that gives the beginning AO index for each atom
c     in the total molecule. 
c     This subroutine makes some assumptions about basis set storing
c     that do not use "legally available" data.
c     Only unilayer runs make sense.
c
      iao=1
      do iat=1,natfmo
        iz=int(fmozan(iat)+0.1D+00)
        ibas=izbas(iat)
        ilay=1
        ist=libish(iz,ibas,ilay)
        ifi=ist+libnsh(iz,ibas,ilay)-1
        nao=0 
        do ii=ist,ifi
          nao=nao+lmax(ii)-lmin(ii)+1
        enddo
        loctat(iat)=iao
        iao=iao+nao 
      enddo
c     write(6,*) 'wwwglol',(loctat(i),i=1,natfmo)
      return
      end
C*MODULE fmolib  *DECK fmomul
      SUBROUTINE fmomul
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      logical GOPARR,DSKWRK,MASWRK
      PARAMETER (MXATM=2000)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /FMCOM / X(1)
      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
c
c     Change multiplicity for FMO-CI runs (monomers+dimers).
c
      mul0=mul
      if(ifmostp.eq.2) then
        mul=ixftch(x(lmulfg),icurfg)
      else if(ifmostp.eq.4) then
        muli=ixftch(x(lmulfg),icurfg)
        mulj=ixftch(x(lmulfg),jcurfg)
        mul=max(muli,mulj) 
c       Such choice is intended to work in the following cases:
c       Singlet and non-singlet -> choose non-singlet 
c       non-singlet and non-singlet -> choose non-singlet 
c       The latter case should mostly be limited to muli=mulj.
      endif
      if(mul.ne.mul0.and.maswrk) write(iw,9000) mul0,mul
      return
 9000 format(/1x,'Warning: multiplicity was changed from',I2,' to',I2/)
      end
C*MODULE fmolib  *DECK pullfrg
      SUBROUTINE pullfrg(ifg0,jndat,ilayfrg,ifgdat,jzbas,indat,iabdfg,
     *                   jabdfg,fmozan,fmoc,izbas,ichfg,mulfg,frgnam,
     *                   layfrg,scffrg,mconfg,molfrg,nprfrg)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      character*8 tstring
      logical GOPARR,DSKWRK,MASWRK
      parameter(MXATM=2000,ncaps=1,maxabd0=4-1,one=1.0D+00,zero=0.0D+00)
      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 /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      dimension jndat(*),ilayfrg(*),ifgdat(*),jzbas(*),indat(*),
     *          iabdfg(*),jabdfg(*),fmozan(*),fmoc(3,*),izbas(*),
     *          ichfg(*),mulfg(*),frgnam(*),layfrg(*),scffrg(*),
     *          mconfg(*),molfrg(*),nprfrg(*),
     *          zcap(ncaps),ccap(3,maxabd0,ncaps),nhcap(ncaps),
     *          a1(3,3),a2(3,3),a3(3,3),zaxis(3),bond(3)
      EQUIVALENCE (tstring,dstring)
      data zaxis/0,0,1/
      data zcap/6.0D+00/,nhcap/3/,
     *     ccap/-0.9709996570D+00,  1.6818207400D+00,-0.6866004421D+00,
     *          -0.9709996570D+00, -1.6818207400D+00,-0.6866004421D+00,
     *           1.9419993140D+00,  0.0000000000D+00,-0.6866004421D+00/
c
c     "Pull out" fragment ifg0, that is, create its "free monomer" data set
c     for RUNTYP=FMO0. 
c
c     capsc stores the library for ncaps prestored sets. It is assumed that
c     all capping atoms are hydrogens, and the atoms to which they are capped
c     are taken from the given molecule. For example, a molecule may have
c     a fractioned bond Ca|---Cb, so for the fragment containing Ca,
c     Cb will be added, and then 3 capping hydrogens to have Ca---CbH3.  
c     At present only CH3 is stored, for R(C-H)=1.09 A. 
c     The central atom should be at the origin. The omitted hydrogen (from CH4) 
c     was put on the z-axis in the positive direction.
c     ifgdat is used to store mapping global indices -> local
c
c     The basis set for the caps is that of the dangling atom.
c     In rare cases the basis set of H may have not be required for normal FMO
c     but will be needed in FMO0. 
c
c     First add the fragment itself.
c
      PI=ACOS(-ONE)
      nat=0
      nfg0=nfg
      nfg=1
      do iat=1,natfmo
        if(indat(iat).eq.ifg0) then
          nat=nat+1
          if(nat.gt.MXATM) call abrt
          zan(nat)=fmozan(iat) 
          jzbas(nat)=izbas(iat)
          call dcopy(3,fmoc(1,iat),1,c(1,nat),1)
          jndat(nat)=nfg
          ifgdat(iat)=nat
        else
          ifgdat(iat)=0
        endif
      enddo
      nat0=nat
c
c     Process fragment-specific options.
c
      ichfg(1)=ichfg(ifg0)
      mulfg(1)=mulfg(ifg0)
      frgnam(1)=frgnam(ifg0)
      ilayfrg(1)=layfrg(ifg0)
      scffrg(1)=scffrg(ifg0)
      mconfg(1)=mconfg(ifg0)
      molfrg(1)=molfrg(ifg0)
      nprfrg(1)=nprfrg(ifg0)
c 
c     write(6,*) 'Added own',nat,nfg,nbdfg
      nbdfgi=0
c     Now add methyl caps. This consists of:
c     adding dangling atoms. Put each into its own fragment.
      do ibdg=1,nbdfg
        ia=abs(iabdfg(ibdg))
        ja=abs(jabdfg(ibdg))
        ifg=indat(ia)
        jfg=indat(ja)
        if(ifg.eq.ifg0) then
          iat=ia
          jat=ja
        else if(jfg.eq.ifg0) then
          iat=ja
          jat=ia
        else
          iat=0
          jat=0
        endif
        if(iat.ne.0) then
c         add the caps. 
c
c         Search the library
          icap=0
          zj=fmozan(jat)
          do i=1,ncaps
            if(abs(zcap(i)-zj).lt.1.0D-06) then
              icap=i
              goto 100
            endif 
          enddo 
          write(6,*) 'Atom not found in the cap library',zj
          call abrt
  100     continue 
c         First add the broken bond information (overwriting the arrays!).
          nat=nat+1
          ialoc=ifgdat(iat)
          jaloc=nat
c         ifgdat should be set for all atoms in ifg0, and iat is guaranteed to 
c         be be in ifg0, so the condition below should never be met. 
          if(ialoc.eq.0) call abrt
          nbdfgi=nbdfgi+1
          if(ifg.eq.ifg0) then
            iabdfg(nbdfgi)=-ialoc
            jabdfg(nbdfgi)= jaloc
          else
            iabdfg(nbdfgi)=-jaloc
            jabdfg(nbdfgi)= ialoc
          endif
          if(nat.gt.MXATM) call abrt
          zan(nat)=fmozan(jat) 
          jzbas(nat)=izbas(jat)
          call dcopy(3,fmoc(1,jat),1,c(1,nat),1)
          nfg=nfg+1
          jndat(nat)=nfg
c
c         Frgament specific options.
c         Make uncharged singlets.
c
          ichfg(nfg)=0
          mulfg(nfg)=1
          WRITE(UNIT=tstring,FMT='(A3,I5.5)') 'cap',ibdg
          frgnam(nfg)=dstring
c         The cap's layer is determined by the layer of the dangling atom. 
          if(ifg.eq.ifg0) then
            ilayfrg(nfg)=layfrg(jfg)
          else
            ilayfrg(nfg)=layfrg(ifg)
          endif
          scffrg(nfg)=scftyp
          mconfg(nfg)=mconfg(1)
          molfrg(nfg)=0
          nprfrg(nfg)=0
c
          call vsub(fmoc(1,jat),1,fmoc(1,iat),1,bond,1,3)
c         write(6,*) 'Added neighbour',nat,nfg,nbdfg
          call vecrot(bond,zaxis,a1)
c
c         Find the closest neighbour of iat (excluding jat). 
c         Search only among ifg0's own atoms (excluding all caps to insure 
c         maximum invariance about input permutations and such).
c         Save time by computing square of the distance.
c
          kat=0
          rk=1.0D+30
          do i=1,nat0
            if(i.ne.ialoc.and.i.ne.jaloc) then
              ri=(c(1,i)-c(1,ialoc))**2 + (c(2,i)-c(2,ialoc))**2 +
     *           (c(3,i)-c(3,ialoc))**2
              if(ri.lt.rk) then
                rk=ri
                kat=i
              endif
            endif
          enddo
c         write(6,*) 'Backbone partner was ',kat
c         if no atom was found (kat=0) it should only mean that there was only
c         one atom in the fragment, in which case the angle does not matter
c         if there is only one cap. In any case recover by setting the angle 
c         to zero. This case of just one atom per frg bonded to others is 
c         probably not practically possible (and using zero angle is not very
c         bad either). 
          if(nat0.eq.1) then
            ANGLE=zero 
            signa=one
          else
            if(kat.eq.0) call abrt
c
c           Add the first capping hydrogen in order to obtain the dihedral
c           angle. Note that nat is not yet modified.
c
            lat=nat+1 
            call MRARBR(a1,3,3,3,ccap(1,1,icap),3,1,c(1,lat),3)
            call daxpy(3,one,fmoc(1,jat),1,c(1,lat),1)
            call DIHED(C,lat,jaloc,ialoc,kat,ANGLE,SIGNA)
c           write(6,*) 'Dihedrals',lat,jaloc,ialoc,kat,SIGNA*ANGLE
          endif
c
c         Build the rotation matrix along the bond, with the angle
c         such that the dihedral angle becomes PI.
c         bond does have to be normalised to call MATNOM.
c         The order of a2*a1 is important! 
c         First rotate R by a1 (rotate z-axis for the caps). 
c         Then rotate a1*R by a2 (rotate caps along the bond).
c
          call dscal(3,one/sqrt(ddot(3,bond,1,bond,1)),bond,1)
          call MATNOM(pi-SIGNA*ANGLE,bond,a2)
          call MRARBR(a2,3,3,3,a1,3,3,a3,3)
c
c         Rotate the capping hydrogens and shift to fmoc(1,jat). 
c
          do i=1,nhcap(icap)
            nat=nat+1
            if(nat.gt.MXATM) call abrt
            zan(nat)=one
            jzbas(nat)=izbas(jat)
            call MRARBR(a3,3,3,3,ccap(1,i,icap),3,1,c(1,nat),3)
            call daxpy(3,one,fmoc(1,jat),1,c(1,nat),1)
            jndat(nat)=nfg
c           write(6,*) 'Added neighbour caps',nat,nfg
          enddo
c         if(kat.ne.0) call DIHED(C,lat,jaloc,ialoc,kat,ANGLE,SIGNA)
c         write(6,*) 'Rotated Dihedral',SIGNA*ANGLE
          if(maswrk) write(iw,9010) ifg0,ibdg,icap,nhcap(icap)+1,kat,
     *                              SIGNA*ANGLE
        endif
      enddo
      nbdfg=nbdfgi
c
c     Rewrite the coordinates+charges.
c
      if(nat.gt.natfmo.or.nfg.gt.nfg0) then
        write(iw,9200) nat,natfmo,nfg,nfg0
        call abrt
      endif
      call dcopy(nat,zan,1,fmozan,1)
      call dcopy(nat*3,c,1,fmoc,1)
      call icopy(nat,jndat,1,indat,1)
      call icopy(nat,jzbas,1,izbas,1)
      call icopy(nfg,ilayfrg,1,layfrg,1)
      natfmo=nat
c     write(6,6666) (zan(kk),(c(jj,kk),jj=1,3),kk=1,nat)
c6666 format(100(F10.5,3F12.5,/))
c     write(6,*) 'bonds',(iabdfg(i),jabdfg(i),i=1,nbdfg)
c
      return
 9010 format(1x,'Free monomer',I5,', ibdg=',I5,', icap=',I2,
     *          ', natc=',I2,', iatdih=',I4,', angle=',F10.6/)
 9200 format(/1x,'Oh no! Your tiny molecule generated a larger capped ',
     *           'system than the original',/1x,'one. The solution is ',
     *           'to add some dummy non-bonded fragment(s).'
     *       /1x,'We suggest one or two H2 coinciding with some real ',
     *           'atoms. Such dummy will have',
     *       /1x,'-zero- effect upon the results but will allocate ',
     *           'enough memory.',
     *       /1x,'You will have to manually remove dummy fragments ',
     *           'from the sparse cube file',
     *       /1x,'though, if you asked for one.',4I5/)
      end
C*MODULE fmolib  *DECK nzCOPY
      SUBROUTINE  nzCOPY(N,DX,INCX,DY,INCY)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION DX(n,INCX),DY(n,INCY)
C
C     COPIES non-zero elements of a VECTOR.
C           DY(I) <== DX(I)
c
      do i=1,n
        if(dx(i,1).ne.0) dy(i,1)=dx(i,1)
      enddo
      return
      end
C*MODULE fmolib  *DECK forcedir 
      logical function forcedir(memory,l1dir,l1,out)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      logical out
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PCKLAB/ LABSIZ
c
c     This function compares the maximum basis set size appropriate for
c     conventional SCF with the current size l1 and returns true if
c     direct SCF should be enfiorced due to insufficient memory.
c
c     The amount is computed based on two assumptions:
c     1.  l1dir (typically, 150-200) basis functions produce 1d+8 integrals,
c     2.  the change to this amount for a different L1 (basis set size) is
c         given by (L1/LDIR)**powerint 
c     Tests indicate powerint=1/3.6 ... 1/3.8, however, we take a bit 
c     aggressive stance here and use 1/3.5. It is still better to dump a bit of
c     integrals on disk rather than recompute them at each SCF iteration. 
c     In addition, the larger the fragments, the lower powerint can get. 
c
c     Memory is the amount per group (not per node), assuming AOINT=DIST 
c     (will not work for sequential MP2?), although this can be
c     adjusted by an appropriate calling argument.
c
c     Input: memory,l1dir,l1,out.
c
      if(l1dir.eq.0.or.memory.eq.0) then
        forcedir=.false.
        return
      else
        faclab=(2+labsiz)/2.0D+00
        powerint=1.0D+00/3.5D+00
c       memory can be negative if 
c       a) DIRSCF=.t. in $SCF and b) NINTIC<0 in $INTGRL.
c       This is a stupid combination, so recover by taking abs of memory. 
c       (For DIRSCF=.f. NINTIC is converted to a positive number earlier).
        maxAOdir=int(l1dir * (abs(memory)/(1.0d+08*faclab))**powerint)
        forcedir=l1.gt.maxAOdir
        if(out) write(iw,9000) maxAOdir,l1,forcedir
      endif
      return
 9000 format(1x,'Max L1 incore=',I6,' current l1=',I6,' force dirscf=',
     *          L2) 
      end
c
C*MODULE fmolib  *DECK dmexch
      SUBROUTINE dmexch(ilay,nfge,jobgrp,layfrg,numfrg,natfrg,scffrg,
     *                  LDAF,iodfmo,maxl30,emon,glocon,nconv,irec0,dosap
     *                 ,dospc,orbxch,enexch,atonce,mastid,popmul,popmat,
     *                  ifgbuf,nfgbuf)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical GOPARR,DSKWRK,MASWRK,dosap,dospc,orbxch,enexch,atonce,
     *        isgddi,parout
      Integer ddi_world,ddi_group
      Parameter(ddi_world=0,ddi_group=1)
      COMMON /FMCOM / X(1)
      COMMON /GDDI/   ISCOPE,NGROUPS,MYGROUP,MEGLOB,npglob,nnglob,JBTYP,
     *                ISGDDI,PAROUT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /TIMING/ CPU,WALL
      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
      dimension jobgrp(*),layfrg(*),numfrg(*),scffrg(*),iodfmo(*),
     *          maxl30(*),emon(*),glocon(*),mastid(0:*),popmul(maxl1,*),
     *        natfrg(*),popmat(maxnat,*),ifgbuf(2,nfg,0:*),nfgbuf(2,0:*)
      data RMC/8HMCSCF   /
c
c       now exchange the initial/current density matrices.
c       the number of converged fragments (nconv) is summed over nodes as
c       part of jobgrp to save one GSUM. 
cnb     skip this subroutine if not running in parallel.
c
c       If massively parallel computing is seriously considered,
c       to increase efficiency one can read and broadcast at once all fragments
c       a group did.
c     parstat: GlobalBcast
c
      CALL TSECND(TIMe0)
      wall0=wall 
      call timit(1)
      if(MASWRK) write(6,*) 'Entering dmexch.'
c     call gddi_scope(ddi_masters)
      jobgrp(nfg+1)=nconv 
      if(.not.maswrk) then
cnb      these are probably not needed
         call viclr(jobgrp,1,nfg+1) 
c        It is wrong to exchange energy for restarts since it will be
c        double counted (duplicated on ALL nodes).
c        if(nconv.ge.0.or.irststp.lt.2) call vclr(emon,1,nfge) 
         if(nconv.ge.0.or.irststp.lt.2) call vclr(emon,1,nfg*2) 
      endif
      mygrp=mygroup
      call gddi_scope(ddi_world)
      ndata=0
c     if(maswrk) then
        CALL DDI_GSUMI(2414,jobgrp,nfg+1)
        nconv=jobgrp(nfg+1)
        if(nconv.ge.0.or.irststp.lt.2) call ddi_gsumf(2422,emon,nfge)
        if(nconv.ge.0) call ddi_gsumf(2422,glocon,ngroups)
        if(maswrk.and.iand(nprfmo,3).le.0) 
     *    write(6,*) 'jobs:',(jobgrp(i),i=1,nfg)
        if(atonce) then
          call viclr(nfgbuf,1,2*ngroups)
          do 100 ifg=1,nfg
            if(layfrg(ifg).lt.ilay) goto 100
            mefg=jobgrp(ifg)
            l1=iand(numfrg(ifg),65535)
            m2=(l1*l1+l1)/2+l1*l1+l1*3 
c           do an upper estimate of m2
            nfgbuf(2,mefg)=nfgbuf(2,mefg)+m2
  100     continue 
          maxb=0
          do i=0,ngroups-1
            maxb=max(maxb,nfgbuf(2,i))
          enddo
c         lbuf1 is used to accumulate data for mygroup 
c         lbuf2 is used for other groups 
          CALL VALFM(LOADFM)
          lbuf1=LOADFM+1
          lbuf2=lbuf1+maxb
          last=lbuf2+maxb
          NEED = LAST- LOADFM -1
          CALL GETFM(NEED)
          ibuf=lbuf1
          call viclr(nfgbuf,1,2*ngroups)
        else
          ibuf=lfmoda
        endif
        do 110 ifg=1,nfg
          if(layfrg(ifg).lt.ilay) goto 110
          l1=iand(numfrg(ifg),65535)
          l2=(l1*l1+l1)/2
          mefg=jobgrp(ifg)
          m2=l2
          if(orbxch) m2=l1*l1
          if(scffrg(ifg).eq.rmc) m2=l2+l1*l1
          if(enexch) m2=m2+l1
          imxl30=maxl30(ifg)
c         The master of the group that did that fragment reads the data in.
          ididit=mastid(mefg)
          nati=natfrg(ifg)
          l1pop=0
          if(dosap) l1pop=l1
          if(dospc) l1pop=l1pop+nati
          m2a=m2+l1pop-nati
          m2t=m2+l1pop
          if(me.eq.ididit) then
            call rareads(LDAF,iodfmo,x(ibuf),m2,ifg+irec0,0)
            if(dosap) call dcopy(l1,popmul(1,ifg),1,x(ibuf+m2),1)
            if(dospc) call dcopy(nati,popmat(1,ifg),1,x(ibuf+m2a),1)
            if(atonce) ibuf=ibuf+m2t
          endif
          ndata=ndata+m2t
          if(atonce) then
c           accumulate all data in lbuf along with indexing 
c           note that indexing must be prepared on all nodes.
            ind=nfgbuf(1,mefg)
            ind=ind+1
            nfgbuf(1,mefg)=ind 
            nfgbuf(2,mefg)=nfgbuf(2,mefg)+m2t
            ifgbuf(1,ind,mefg)=ifg
            ifgbuf(2,ind,mefg)=m2
          else
            CALL DDI_BCAST(2416,'F',x(ibuf),m2t,ididit)
c           it may be a bit of a philosphical decision: whether to enforce
c           having the same density for each node here. At present group
c           members do not write what their master broadcast because they
c           should have the same thing anyway (so we think).
c           if(mygrp.ne.mefg) then
            if(me.ne.ididit) then
              CALL rawrites(LDAF,iodfmo,x(ibuf),imxl30,m2,ifg+irec0,0)
              if(dosap) call dcopy(l1,x(ibuf+m2),1,popmul(1,ifg),1)
              if(dospc) call dcopy(nati,x(ibuf+m2a),1,popmat(1,ifg),1)
            endif
          endif
c         if(dosap.and.maswrk) then
c           write(6,*) ifg,'orbital Mullikens',l1
c           call prsq(popmul(1,ifg),l1,1,1)
c         endif
c         if(dospc.and.maswrk) then
c           write(6,*) ifg,'atomic Mullikens',nat0
c           call prsq(popmat(1,ifg),nat0,1,1)
c         endif
  110   continue 
        if(atonce) then
c       now broadcast the accumulated data
          do i=0,ngroups-1
            m2t=nfgbuf(2,i) 
            if(m2t.ne.0) then
            ididit=mastid(i)
c           sending side provides mygroup data in lbuf1, receiving side
c           stores to lbuf2, in order not to overwrite its own data.
            ibuf=lbuf2
            if(mygrp.eq.i) ibuf=lbuf1
            CALL DDI_BCAST(2416,'F',x(ibuf),m2t,ididit)
c           if(me.ne.ididit) then
            if(mygrp.ne.i) then
              nn=nfgbuf(1,i)
              do j=1,nn
                ifg=ifgbuf(1,j,i)
                m2= ifgbuf(2,j,i)
                imxl30=maxl30(ifg)
                l1=iand(numfrg(ifg),65535)
                nati=natfrg(ifg)
                l1pop=0
                if(dosap) l1pop=l1
                if(dospc) l1pop=l1pop+nati
                m2a=m2+l1pop-nati
                CALL rawrites(LDAF,iodfmo,x(ibuf),imxl30,m2,ifg+irec0,0)
                if(dosap) call dcopy(l1,x(ibuf+m2),1,popmul(1,ifg),1)
                if(dospc) call dcopy(nati,x(ibuf+m2a),1,popmat(1,ifg),1)
                ibuf=ibuf+m2+l1pop
              enddo
            endif
            else
              if(maswrk) write(6,*) 'Nothing done by group',i
            endif
          enddo
          call retfm(need) 
        endif
c     endif
      call gddi_scope(ddi_group)
c     CALL DDI_BCAST(2417,'I',nconv,1,master)
      call timit(1)
      CALL TSECND(TIMe1)
      if(maswrk) write(iw,9600) ndata,wall-wall0
      return
 9600 format(/1x,'Data exchange of',I10,' words took',F8.1,' s.')
      end
c
C*MODULE fmolib  *DECK fmopre
      subroutine fmopre(ilay,nstep,ngrfmo,mannod,mastid)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (maxpst=10)
      dimension ngrfmo(maxpst,*),mannod(*),mastid(*)
c
      call timit(1)
c
c     Find the offset: mannod is stored as an array of vectors with variable
c     size. Each vector contains one set of group sizes.
c     This is repeated for each layer.
c
      ioff=1
      do i=1,ilay-1
        do j=1,maxpst
          ioff=ioff+ngrfmo(j,i)
        enddo
      enddo
      do i=1,nstep-1
        ioff=ioff+ngrfmo(i,ilay)
      enddo
      call gddi_init(ngrfmo(nstep,ilay),mannod(ioff))
      call gddi_mastid(mastid)
c
      call timit(1)
c
      return
      end
