c 17 Jan 06 - DGF - FMOCHECK - alter FMO/PCM test
C 21 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 31 Jan 05 - DGF,KK - add property module for the FMO method
c
C*MODULE fmoprp  *DECK mongrid
      SUBROUTINE mongrid(NGm,ORIGM,UXM,UYM,UZM,vdwr,pad)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter (MXATM=2000)
      logical GOPARR,DSKWRK,MASWRK
      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
      dimension NGm(3),ORIGM(3),UXM(3),UYM(3),UZM(3),vdwr(*),Cmin(3),
     *          Cmax(3)
c
c     Generate grid automatically. 
c     At present a rectangular grid is constructed surrounding the molecule
c     and leaving an extra margin of pad vdW radii from all sides.
c
      if(NGm(1).le.1.or.NGm(2).le.1.or.NGm(3).le.1) call abrt
c
      call dacopy(3, 1.0D+20,Cmin,1)
      call dacopy(3,-1.0D+20,Cmax,1)
      do i=1,nat
        rr=vdwr(ian(i))*pad 
        do j=1,3
          Cmin(j)=min(Cmin(j),C(j,i)-rr)
          Cmax(j)=max(Cmax(j),C(j,i)+rr)
        enddo
      enddo
      call vclr(UXM,1,3) 
      call vclr(UYM,1,3) 
      call vclr(UZM,1,3) 
      call dcopy(3,Cmin,1,ORIGM,1)
      UXM(1)=(Cmax(1)-Cmin(1))/(NGm(1)-1)
      UYM(2)=(Cmax(2)-Cmin(2))/(NGm(2)-1)
      UZM(3)=(Cmax(3)-Cmin(3))/(NGm(3)-1)
      if(maswrk) write(iw,9000) Cmin(1),Cmin(2),Cmin(3),
     *                          Cmax(1),Cmax(2),Cmax(3)
      return
 9000 format(/1x,'Grid is from',3F10.4,' to ',3F10.4)
      end
c
C*MODULE fmoprp  *DECK grdbox
      SUBROUTINE grdbox(ixmin,ixmax,iymin,iymax,izmin,izmax,vdwr,pad)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter (MXATM=2000)
      logical GOPARR,DSKWRK,MASWRK
      COMMON /GRDPAR/ ORIGIN(3),XVEC(3),YVEC(3),ZVEC(3),UX(3),UY(3),
     *                UZ(3),GRDSIZ,NGRID,IGUNIT,NXG,NYG,NZG,modgrd
      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
      dimension vdwr(*),Cmin(3),Cmax(3)
c
c     Generate a box surrounding the current n-mer for the purpose of
c     computing electron density on grid. 
c     At present a rectangular grid is constructed surrounding the molecule
c     and leaving an extra margin of pad vdW radii from all sides.
c
      uxl=sqrt(ddot(3,ux,1,ux,1))
      uyl=sqrt(ddot(3,uy,1,uy,1))
      uzl=sqrt(ddot(3,uz,1,uz,1))
      uxo=ddot(3,origin,1,ux,1)/uxl 
      uyo=ddot(3,origin,1,uy,1)/uyl
      uzo=ddot(3,origin,1,uz,1)/uzl 
c
      call dacopy(3, 1.0D+20,Cmin,1)
      call dacopy(3,-1.0D+20,Cmax,1)
c     Find the smallest and the largest projections of the molecule upon the 
c     grid axes. The extra space due to vdW radii is taken slightly over,
c     by surrounding the projection unto axes rather than formally drawing a 
c     sphere and finding the maximum and minimum projection of a sphere.
      do i=1,nat
        rr=vdwr(ian(i))*pad
        uxp=(C(1,i)*ux(1)+C(2,i)*ux(2)+C(3,i)*ux(3))/uxl
        uyp=(C(1,i)*uy(1)+C(2,i)*uy(2)+C(3,i)*uy(3))/uyl
        uzp=(C(1,i)*uz(1)+C(2,i)*uz(2)+C(3,i)*uz(3))/uzl
        Cmin(1)=min(Cmin(1),uxp-rr)
        Cmax(1)=max(Cmax(1),uxp+rr)
        Cmin(2)=min(Cmin(2),uyp-rr)
        Cmax(2)=max(Cmax(2),uyp+rr)
        Cmin(3)=min(Cmin(3),uzp-rr)
        Cmax(3)=max(Cmax(3),uzp+rr)
      enddo
      ixmin=max(int((Cmin(1)-uxo)/uxl)+1,1)
      ixmax=min(int((Cmax(1)-uxo)/uxl)+2,NXG)
      iymin=max(int((Cmin(2)-uyo)/uyl)+1,1)
      iymax=min(int((Cmax(2)-uyo)/uyl)+2,NYG)
      izmin=max(int((Cmin(3)-uzo)/uzl)+1,1)
      izmax=min(int((Cmax(3)-uzo)/uzl)+2,NZG)
c
c     guard against "negative boxes", that arise if some fragments have
c     nothing in the total molecular box (with manual box given in $GRID).
c     A way to recover is to use the (1,0) window in that dimension, which
c     gives 0 points (imax-imin+1). Note that memory allocation will be correct
c     based on 0 points in a grid window and everything else is supposed to
c     function properly.
c
      iempty=0
      if(ixmin.le.0.or.ixmax.le.0.or.ixmin.gt.ixmax) then
        ixmin=1
        ixmax=0
        iempty=1
      endif
      if(iymin.le.0.or.iymax.le.0.or.iymin.gt.iymax) then
        iymin=1
        iymax=0
        iempty=2
      endif
      if(izmin.le.0.or.izmax.le.0.or.izmin.gt.izmax) then
        izmin=1
        izmax=0
        iempty=3
      endif
      if(maswrk) then
        if(iempty.eq.0) then
          write(iw,9000) ixmin,ixmax,iymin,iymax,izmin,izmax 
        else
          write(iw,9010) ixmin,ixmax,iymin,iymax,izmin,izmax 
        endif
      endif
      return
 9000 format(1x,'Reduced grid to',6I5)
 9010 format(1x,'Reduced grid to',6I5,' : empty box!')
      end
c
C*MODULE fmoprp  *DECK totgrid
      SUBROUTINE totgrid(fmozan,fmoc,vdwr,pad)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      logical GOPARR,DSKWRK,MASWRK
      COMMON /GRDPAR/ ORIGIN(3),XVEC(3),YVEC(3),ZVEC(3),UX(3),UY(3),
     *                UZ(3),GRDSIZ,NGRID,IGUNIT,NXG,NYG,NZG,modgrd
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      dimension fmozan(*),fmoc(3,*),vdwr(*),Cmin(3),Cmax(3)
c
c     Generate grid automatically.
c     At present a rectangular grid is constructed surrounding the molecule
c     and leaving an extra margin of pad vdW radii from all sides.
c
      call dacopy(3, 1.0D+20,Cmin,1)
      call dacopy(3,-1.0D+20,Cmax,1)
      do i=1,natfmo
        iz=int(fmozan(i)+0.1D+00)
        rr=vdwr(iz)*pad
        do j=1,3
          Cmin(j)=min(Cmin(j),fmoc(j,i)-rr)
          Cmax(j)=max(Cmax(j),fmoC(j,i)+rr)
        enddo
      enddo
      call vclr(UX,1,3)
      call vclr(UY,1,3)
      call vclr(UZ,1,3)
      call dcopy(3,Cmin,1,ORIGIN,1)
c     if(GRDSIZ.ne.0.0D+00) then
        UX(1)=GRDSIZ
        UY(2)=GRDSIZ
        UZ(3)=GRDSIZ
        NXG=int((Cmax(1)-Cmin(1))/GRDSIZ+0.5D+00)
        NYG=int((Cmax(2)-Cmin(2))/GRDSIZ+0.5D+00)
        NZG=int((Cmax(3)-Cmin(3))/GRDSIZ+0.5D+00)
        if(maswrk) write(iw,9010) GRDSIZ,NXG,NYG,NZG
c     else
c       if(NXG.le.1.or.NYG.le.1.or.NZG.le.1) call abrt
c       UX(1)=(Cmax(1)-Cmin(1))/(NXG-1)
c       UY(2)=(Cmax(2)-Cmin(2))/(NYG-1)
c       UZ(3)=(Cmax(3)-Cmin(3))/(NZG-1)
c     endif
      if(maswrk) write(iw,9000) Cmin(1),Cmin(2),Cmin(3),
     *                          Cmax(1),Cmax(2),Cmax(3)
      return
 9000 format(/1x,'Grid is from',3F16.8,
     *       /10x,        'to ',3F16.8,' (bohr)',/)
 9010 format(/1x,'SIZE=',F10.4,' bohr, ',I5,' by ',I5,' by ',I5,
     *           ' grid points.') 
      end
c
C*MODULE fmoprp  *DECK makecube
      SUBROUTINE makecube(ifg,ibody,natf,fmozan,fmoc,dgrid,
     *                    NXG,NYG,NZG,ORIGIN,UX,UY,UZ) 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      parameter (tiny=1.0d-06)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
c     Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      DIMENSION fmozan(*),fmoc(3,*),dgrid(NZG,NYG,NXG),
     *          ORIGIN(3),UX(3),UY(3),UZ(3)
c
c     Create Gaussian-style cube file with numeric properties on a grid.
c
      if(ifg.eq.0) then
        write(ip,9000) ibody
      else
c       write(ip,9005) min(nbody,1),ifg
        write(ip,9005) ibody,ifg
      endif
c
c     the sign of iflag means: +1 density, -1 MOs 
c
      iflag=natf
c
      write(ip,9010) iflag,(ORIGIN(i),i=1,3)
      write(ip,9010) NXG,(UX(i),i=1,3)
      write(ip,9010) NYG,(UY(i),i=1,3)
      write(ip,9010) NZG,(UZ(i),i=1,3)
      do i=1,natf
        write(ip,9010) int(fmozan(i)+tiny),fmozan(i),(fmoc(j,i),j=1,3)
      enddo
      do ix=1,NXG
        do iy=1,NYG
          write(ip,9020) (dgrid(iz,iy,ix),iz=1,NZG)
        enddo 
      enddo 
      write(ip,9100)
      return
 9000 format(1x,'$CUBE',I1,
     *      /1x,'DENSITY: full system, created by GAMESS (FMO).',/)
 9005 format(1x,'$CUBE',I1,
     *      /1x,'DENSITY: fragment',I5,', created by GAMESS (FMO).',/)
 9010 format(i5,4f12.6)
 9020 format(6e13.5)
 9100 format(1x,'$END')
      end
c
C*MODULE fmoprp  *DECK makescube
      SUBROUTINE makescube(mode,ifgfmo0,ibody,layfrg,indat,fmozan,fmoc,
     *                     dgrid,noffg,indgrd,maxg)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      CHARACTER*30 fmtstr
      parameter (tiny=1.0d-06)
      COMMON /GRDPAR/ ORIGIN(3),XVEC(3),YVEC(3),ZVEC(3),UX(3),UY(3),
     *                UZ(3),GRDSIZ,NGRID,IGUNIT,NXG,NYG,NZG,modgrd
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      DIMENSION layfrg(*),indat(*),fmozan(*),fmoc(3,*),dgrid(*),noffg(*)
     *         ,indgrd(6,*),maxg(3)
c
c     Create sparse cube file with numeric properties on a grid.
c     mode=0 only header
c     mode=1 header+data
c
c     the sign of iflag means: +1 density, -1 MOs 
c
      iflag=natfmo

      if(ifgfmo0.eq.0) then
        write(ip,9000) ibody
      else
        if(mode.eq.0) then
          igrp=0
          nfingers=1
        else
          igrp=ifgfmo0
          fgrp=igrp
          nfingers=int(log10(fgrp))+1
        endif
        if(ifgfmo0.eq.1.or.mode.ne.0) then 
          WRITE(UNIT=fmtstr,FMT='(12H(1x,''$SCM'',I,I1,1H))') nfingers
          WRITE(ip,FMT=fmtstr) igrp
        endif 
        if(ifgfmo0.eq.1.and.mode.eq.0) write(ip,9007)
      endif
      if(ifgfmo0.eq.0.or.(ifgfmo0.eq.1.and.mode.eq.0)) then
        write(ip,9010) nfg,iflag,(ORIGIN(i),i=1,3)
        write(ip,9010) maxg(1),NXG,(UX(i),i=1,3)
        write(ip,9010) maxg(2),NYG,(UY(i),i=1,3)
        write(ip,9010) maxg(3),NZG,(UZ(i),i=1,3)
        do i=1,natfmo
          write(ip,9010) indat(i),int(fmozan(i)+tiny),
     *                   fmozan(i),(fmoc(j,i),j=1,3)
        enddo
      endif
c
      if(mode.ne.0) then
      nfg0=nfg
      if(ifgfmo0.ne.0) nfg0=1
      do ifg=1,nfg0
        ioff=noffg(ifg)
        ixmin=indgrd(1,ifg)
        ixmax=indgrd(2,ifg)
        iymin=indgrd(3,ifg)
        iymax=indgrd(4,ifg)
        izmin=indgrd(5,ifg)
        izmax=indgrd(6,ifg)
        NXGi=ixmax-ixmin+1
        NYGi=iymax-iymin+1
        NZGi=izmax-izmin+1
        write(ip,9040) ifg,layfrg(ifg),ixmin,ixmax,iymin,iymax,
     *                 izmin,izmax
        do ix=1,NXGi
          do iy=1,NYGi
            write(ip,9050) (dgrid(ioff+iz-1),iz=1,NZGi)
            ioff=ioff+NZGi
          enddo 
        enddo 
      enddo 
      endif 
      write(ip,9100)
c
      return
 9000 format(1x,'$SCUBE',I1,/1x,'DENSITY: sparse, full system, ',
     *          'created by GAMESS (FMO).',/)
 9007 format(1x,
     *      'DENSITY: sparse, free monomer, created by GAMESS (FMO).',/)
 9010 format(2i6,4f12.6)
 9040 format(8i6)
 9050 format(6e13.5)
 9100 format(1x,'$END')
      end
c
C*MODULE fmoprp  *DECK proindt
      SUBROUTINE proindt(nfg2,nfg3,ilay,loadhf,layfrg,numfrg,molfrg,
     *                   nmtfrg,ndtfrg,iwrk,loadt,net,netc)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      logical nxt,isgddi,parout,goparrm,trimc
      COMMON /GDDI/   ISCOPE,NGROUPS,MYGROUP,MEGLOB,npglob,nnglob,JBTYP,
     *                ISGDDI,PAROUT
      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
      dimension layfrg(*),numfrg(*),molfrg(*),nmtfrg(nfg,2),
     *          ndtfrg(nfg2,2),iwrk(*),loadt(*)
c
c     This subroutine fills multiplicative constant arrays used to compute
c     properties for trimers. This subroutine should be called after global
c     summing the distances so that they are available.
c     Note that at present only masters call this subroutine and it is a 
c     rare example of load balancing done within the master scope.
c     This subroutine ignores reordering of monomers (size load balancing).
c     The results are summed above (where this subroutine is called).
c     Dynamic load balancing is ignored as well because it does not work
c     in the master scope?
c     Note that this subroutine counts trimc assuming correlation. If none
c     is there netc should be zeroed out elsewhere.
c
      call viclr(nmtfrg,1,nfg*2)
      call viclr(ndtfrg,1,nfg2*2)
      rmin=0
      rmax=0
      rrij=0
      rrik=0
      rrjk=0
      ishift=maxl1*3
c     NXT = IBTYP.EQ.1
      NXT = .false.
      NEXT  = -1
      kount = -1
      nprocm=ngroups
      goparrm=nprocm.gt.1
      mem=mygroup 
c     write(6,*) NPROC,goparr,me
c     neither of nproc, goparr, me is set properly for the master scope.
      ijk=0
      do 100 ifg=1,nfg
        kount=kount+1
        if(goparrm) then
          n2skip=(ifg-1)*(ifg-2)
          IF(NXT) THEN
            IF(kount.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
            if(NEXT.ne.kount) then
              ijk=ijk+n2skip/2
              goto 100
            endif
          else
            if(MOD(kount,NPROCm).NE.mem) then
              ijk=ijk+n2skip/2
              goto 100
            endif
          endif
        endif
        ifglay=layfrg(ifg)
        do jfg=1,ifg-1
          jfglay=layfrg(jfg)
          ijfg=(ifg*ifg-3*ifg)/2+jfg+1
          if(resdim.ne.0) rrij=fmodist(ifg,0,0,jfg)
          do kfg=1,jfg-1
            kfglay=layfrg(kfg)
            ijk=ijk+1
            l1ijk=iand(numfrg(ifg),65535)+iand(numfrg(jfg),65535)+
     *            iand(numfrg(kfg),65535)
            loo=l1ijk
            if((ifglay.eq.ilay.or.jfglay.eq.ilay.or.kfglay.eq.ilay).and.
     *        ifglay.ge.ilay.and.jfglay.ge.ilay.and.kfglay.ge.ilay) then
              if(restri(1).ne.0) call fmodist3(ifg,jfg,kfg,rmin,rmax)
              if(((rmin.le.restri(1).or.rmax.le.restri(2)).and.(rmax.le.
     *            restri(3).or.restri(3).eq.0)).or.restri(1).eq.0) then
                if(resdim.ne.0) then
                  rrik=fmodist(ifg,0,0,kfg)
                  rrjk=fmodist(jfg,0,0,kfg)
                endif
                ires=molfrg(ifg)+molfrg(jfg)+molfrg(kfg)
c               no change yet for exclmol
                trimc=(rmax.le.restri(4).or.restri(4).eq.0).and.
     *                (ires.lt.0.or.ires.eq.1) 
                if(rrij.le.resdim.or.resdim.eq.0) then
                  ndtfrg(ijfg,1)=ndtfrg(ijfg,1)+1
                  nmtfrg(ifg,1)=nmtfrg(ifg,1)+1
                  nmtfrg(jfg,1)=nmtfrg(jfg,1)+1
                  ires=molfrg(ifg)+molfrg(jfg)
                  if(trimc.and.(rrij.le.rcorsd.or.rcorsd.eq.0).and.
     *               ires.ne.0.and.ires.ne.2) then
                    ndtfrg(ijfg,2)=ndtfrg(ijfg,2)+1
                    nmtfrg(ifg,2)=nmtfrg(ifg,2)+1
                    nmtfrg(jfg,2)=nmtfrg(jfg,2)+1
                  endif
                endif
                if(rrik.le.resdim.or.resdim.eq.0) then
                  ikfg=(ifg*ifg-3*ifg)/2+kfg+1
                  ndtfrg(ikfg,1)=ndtfrg(ikfg,1)+1
                  nmtfrg(ifg,1)=nmtfrg(ifg,1)+1
                  nmtfrg(kfg,1)=nmtfrg(kfg,1)+1
                  ires=molfrg(ifg)+molfrg(kfg)
                  if(trimc.and.(rrik.le.rcorsd.or.rcorsd.eq.0).and.
     *               ires.ne.0.and.ires.ne.2) then
                    ndtfrg(ikfg,2)=ndtfrg(ikfg,2)+1
                    nmtfrg(ifg,2)=nmtfrg(ifg,2)+1
                    nmtfrg(kfg,2)=nmtfrg(kfg,2)+1
                  endif
                endif
                if(rrjk.le.resdim.or.resdim.eq.0) then
                  jkfg=(jfg*jfg-3*jfg)/2+kfg+1
                  ndtfrg(jkfg,1)=ndtfrg(jkfg,1)+1
                  nmtfrg(jfg,1)=nmtfrg(jfg,1)+1
                  nmtfrg(kfg,1)=nmtfrg(kfg,1)+1
                  ires=molfrg(jfg)+molfrg(kfg)
                  if(trimc.and.(rrjk.le.rcorsd.or.rcorsd.eq.0).and.
     *               ires.ne.0.and.ires.ne.2) then
                    ndtfrg(jkfg,2)=ndtfrg(jkfg,2)+1
                    nmtfrg(jfg,2)=nmtfrg(jfg,2)+1
                    nmtfrg(kfg,2)=nmtfrg(kfg,2)+1
                  endif
                endif
                do i=1,2
                  if(i.eq.1.or.trimc) then
                    nmtfrg(ifg,i)=nmtfrg(ifg,i)-1
                    nmtfrg(jfg,i)=nmtfrg(jfg,i)-1
                    nmtfrg(kfg,i)=nmtfrg(kfg,i)-1
                  endif
                enddo
                net=net+1
                loo=loo+ishift
                if(trimc) then
                  netc=netc+1
                  loo=loo+ishift
                endif
              endif 
            endif
            if(loadhf.ne.0) iwrk(ijk)=loo
          enddo
        enddo
  100 continue
      if(goparrm.and.nxt) CALL DDI_DLBRESET
      if(loadhf.ne.0) then
        call ddi_gsumi(2421,iwrk,nfg3)
        call indsort(nfg3,iwrk,loadt)
c       write(6,*) 'nloadt',(loadt(i),i=1,nfg3)
      endif
      RETURN
      END
c
C*MODULE fmoprp  *DECK fmodeg
      subroutine fmodeg(ida,fmode,iaglob)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension fmode(3,*),iaglob(*)
      Parameter (MXATM=2000)
      common /GRAD  / DE(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
c
c     gather FMO gradient contributions coming from the analytic gradient of a
c     monomer/dimer.
c
c     write(6,*) 'Grad accs',ida
      if(ida.eq.0) return
c     CALL EGOUT(DE,NAT)
      da=ida
      do i=1,nat
        ig=iaglob(i)
        do j=1,3
          fmode(j,ig)=fmode(j,ig)+da*de(j,i)
        enddo
      enddo
      return
      end
c
C*MODULE fmoprp  *DECK monor
      subroutine monor(ifg,indat,fmozan,fmoc,iabdfg,jabdfg,rij,rwrk,vdwr
     *                )
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension indat(*),fmozan(*),fmoc(3,*),iabdfg(*),
     *          jabdfg(*),vdwr(*),rij(*),rwrk(*)
      PARAMETER (MXATM=2000,zero=0.0D+00,ten=1.0D+01)
      LOGICAL GOPARR,DSKWRK,MASWRK,nxt
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,convfg,cnvdmp,coroff,IXESP,mxitfg,nguess,
     *                NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c
c
c     fill in interfragment distances relative to van der Walls radii, that is,
c     rIJ=min(i(I),j(J)) {rij/(Ri+Rj)}, where R are van der Waals radii. 
c
      nifg=ifg-1
c     initialise with very large numbers 
      call dacopy(nifg,1.0D+30,rij,1)
c
c     fill in cut bonds zero distances 
c
      do 500 ibdg=1,nbdfg
        ibda=iabdfg(ibdg)
        jbda=jabdfg(ibdg)
        ibfg=indat(abs(ibda)) 
        jbfg=indat(abs(jbda)) 
        if(jbda.lt.0.and.jbfg.eq.ifg) then
          rij(ibfg)=zero
          write(6,*) 'should never come here; check ibdfg and jbdfg'
c         that mean the user has given noncanonical bond cleavage
          call abrt
        endif
        if(jbfg.ge.ifg.or.rij(jbfg).eq.zero) goto 500
        if(ibda.lt.0.and.ibfg.eq.ifg) then
          rij(jbfg)=zero
          goto 500
        endif
c       assume canonical ordering and process ghost atoms
c       Every atom can belong to several fragments. The "main" N-1 fragment
c       is processed in the 200 loop below and here the shattered ghost atoms
c       are treated. This is why coordinates are from the N-1 side (ibda)
c       and it is assigned to the ghost side (jbfg). 
c       if(jbfg.eq.ifg) goto 500
c       Skip ghost atoms within monomer ifg. 
        j=abs(ibda) 
        rj=vdwr(int(fmozan(j)+1.0D-03))
        rr12=1.0D+30
        xj=fmoc(1,j)
        yj=fmoc(2,j)
        zj=fmoc(3,j)
        do i=1,nat
c         ri=vdwr(int(zan(i)+1.0D-03))
          ri=vdwr(IAN(i))
          r12=sqrt((c(1,i)-xj)**2+(c(2,i)-yj)**2+(c(3,i)-zj)**2)/(ri+rj)
          if(r12.lt.rr12) rr12=r12
        enddo
        rij(jbfg)=min(rij(jbfg),rr12)
c       rij(jbfg)=min(rij(jbfg),rr12*UNITS)
  500 continue 
c
      NXT = IBTYP.EQ.1
      NEXT  = -1
      kount = -1
      nratch=3
      if(.not.goparr.or..not.nxt) nratch=1
      nchsiz=(natfmo-1)/(nproc*nratch)+1
c     nratch is a DLB parallelisation ratio. The work is divided so that each 
c     node does NRATCH chunks on the average.  
c     Now the initial value of nratch is hardwired. The value can be
c     reconsidered if a massive number of nodes per group is available
c     (or if nodes of different type are mixed).
c     Static work division is not good because 7 is divided as 3, 3 and 1. 
c
c     divide natfmo into chunks 
c 
      do 200 j1=1,natfmo,nchsiz
        kount=kount+1
        if(goparr) then
          IF(NXT) THEN
            IF(kount.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
            if(NEXT.ne.kount) goto 200
          else
            if(MOD(kount,NPROC).NE.me) goto 200
          endif
        endif
        j2max=min(nchsiz-1,natfmo-j1)
c       note that j2 starts at 0. 
        do 100 j2=0,j2max
          j=j1+j2
          jfg=indat(j)
          if(jfg.ge.ifg.or.rij(jfg).eq.zero) goto 100
          rj=vdwr(int(fmozan(j)+1.0D-03))
          xj=fmoc(1,j)
          yj=fmoc(2,j)
          zj=fmoc(3,j)
          rr12=1.0D+30
          do i=1,nat
c           ri=vdwr(int(zan(i)+1.0D-03))
            ri=vdwr(IAN(i))
          r12=sqrt((c(1,i)-xj)**2+(c(2,i)-yj)**2+(c(3,i)-zj)**2)/(ri+rj)
            if(r12.lt.rr12) rr12=r12
          enddo
          rij(jfg)=min(rij(jfg),rr12) 
c         rij(jfg)=min(rij(jfg),rr12*UNITS) 
  100   continue 
  200 continue 
      if(goparr) then
        IF(NXT) CALL DDI_DLBRESET
c       DDI_GSUMF will not work, we need DDI_GMINF!
c       find maximum over all nodes.
        do 300 my=0,nproc-1
          if(me.eq.my) call dcopy(nifg,rij,1,rwrk,1)
          CALL DDI_BCAST(2422,'F',rwrk,nifg,my)
          if(me.ne.my) then
            do i=1,nifg
              rij(i)=min(rij(i),rwrk(i))
            enddo
          endif
  300   continue  
      endif 
      if(iand(nfmopal,1).ne.0) then
        do i=1,nifg
c         irij=int(rij(i)*ten+5.0D-01)
          irij=int(rij(i)*ten)
          rij(i)=irij/ten
        enddo
      endif
      return
      end
C*MODULE fmoprp  *DECK mulpop
      subroutine mulpop(l1,d,s,popmul)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension d(*),s(*),popmul(*)
      PARAMETER (zero=0.0D+00)
c
c     compute Mulliken orbital populations
c
c     call prtril(d,l1)
c     call prtril(s,l1)
      loop=0
      do i=1,l1
        sum=zero
        do j=1,i-1
          loop=loop+1
          a=d(loop)*s(loop) 
          popmul(j)=popmul(j)+a
          sum=sum+a
        enddo
        loop=loop+1
        popmul(i)=sum+d(loop)*s(loop)
      enddo
      return
      end
C*MODULE fmoprp  *DECK mulpopa
      subroutine mulpopa(popmul,popmat)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000,MXSH=5000,MXGTOT=20000,zero=0.0D+00)
      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 /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      dimension popmul(*),popmat(*)
c
cnb   fix that
c     due to dmexch difficulties, clear out all
c     call vclr(popmat,1,nat)
      call vclr(popmat,1,max(maxnat,nat))
      do ii=1,NSHELL
        iat=katom(ii)
        mini=kmin(ii)
        maxi=kmax(ii)
        loci=kloc(ii)-mini
        sum=zero
        do i=mini,maxi
          sum=sum+popmul(loci+i)
        enddo
        popmat(iat)=popmat(iat)+sum 
      enddo
      return
      end
C*MODULE fmoprp  *DECK addprop
      subroutine addprop(ibody,primulm,primul,idam,idad,idat,iaglob,
     *                   popmat,nfmoelm,fmoelm,fmoelmd,fmoelmt,atmulq,
     *                   atmulqd,atmulqt)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical primulm,primul,GOPARR,DSKWRK,MASWRK
      dimension fmoelm(*),fmoelmd(*),fmoelmt(*),iaglob(*),popmat(*),
     *          atmulq(*),atmulqd(*),atmulqt(*)
      Parameter (MXATM=2000)
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      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 /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /XYZPRP/ XP,YP,ZP,
     *                DMX,DMY,DMZ,
     *                QXX,QYY,QZZ,QXY,QXZ,QYZ,
     *                QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ,
     *                OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                OXZZ,OYZZ,OZZZ,OXYZ,
     *                OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
      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
c
c     Accumulate monomer properties, that is, electrostatic moments,
c     and Mulliken charges.
      dad=idad
      dat=idat
c
c     Compute electrostatic moments. Only the masters do that to avoid double
c      counting in global sums.
c
      if(nfmoelm.ne.0) then
        if(nprint.eq.-5) then
c       Then SCF ignored properties!!
          nprints=nprint
          if(iand(nprfmo,3).le.1) nprint=0
          call ELMOMC
          nprint=nprints
        endif
        if(.not.maswrk) return
        if(iemom.ge.1) then
          fmoelm(1)= fmoelm(1)+DMX
          fmoelm(2)= fmoelm(2)+DMY
          fmoelm(3)= fmoelm(3)+DMZ
        endif
        if(iemom.ge.2) then
          fmoelm(4)= fmoelm(4)+QXX
          fmoelm(5)= fmoelm(5)+QYY
          fmoelm(6)= fmoelm(6)+QZZ
          fmoelm(7)= fmoelm(7)+QXY
          fmoelm(8)= fmoelm(8)+QXZ
          fmoelm(9)= fmoelm(9)+QYZ
        endif
        if(iemom.ge.3) then
          fmoelm(10)=fmoelm(10)+OXXX
          fmoelm(11)=fmoelm(11)+OXXY
          fmoelm(12)=fmoelm(12)+OXXZ
          fmoelm(13)=fmoelm(13)+OXYY
          fmoelm(14)=fmoelm(14)+OYYY
          fmoelm(15)=fmoelm(15)+OYYZ
          fmoelm(16)=fmoelm(16)+OXZZ
          fmoelm(17)=fmoelm(17)+OYZZ
          fmoelm(18)=fmoelm(18)+OZZZ
          fmoelm(19)=fmoelm(19)+OXYZ
        endif
        if(ibody.eq.1.and.nbody.gt.1.or.ibody.eq.2.and.nbody.gt.2)
     *    call daxpy(nfmoelm,dad,fmoelm,1,fmoelmd,1)
        if(ibody.eq.1.and.nbody.gt.2)
     *    call daxpy(nfmoelm,dat,fmoelm,1,fmoelmt,1)
      endif
      if(.not.maswrk) return
c
c     compute Mulliken populations 
c
      if(primul) then
        if(primulm.and.idam.ne.0) then
          do iat=1,nat
            iatg=iaglob(iat)
            atmulq(iatg)=atmulq(iatg)+popmat(iat)
          enddo
        endif
        if(maswrk.and.nprint.eq.-5) then 
          write(iw,9400)
          do iat=1,nat
            write(iw,9410) iat,zan(iat),popmat(iat)
          enddo
        endif
        if(ibody.lt.nbody.and.nbody.gt.1) then 
          do iat=1,nat
            iatg=iaglob(iat)
            atmulqd(iatg)=atmulqd(iatg)+dad*popmat(iat)
          enddo
          if(ibody.eq.1.and.nbody.gt.2) then
            do iat=1,nat
              iatg=iaglob(iat)
              atmulqt(iatg)=atmulqt(iatg)+dat*popmat(iat)
            enddo
          endif
        endif
      endif
      return
 9400 format(/1x,'Mulliken atomic populations',/5x,'I',3x,'Z',4x,'Q')
 9410 format(1x,I5,F5.1,F10.7)
      END
c
C*MODULE fmoprp  *DECK fmocheck
      subroutine fmocheck(nder,mulfg,scffrg,fmoscf,fmoci,fmodft,fmocc,
     *                    mpnfmo)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (zero=0.0D+00)
      logical GOPARR,DSKWRK,MASWRK,domcscf,doci,dodft,docc,domp,
     *        legalscf,legalmul
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEXTFLD
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IRPPCM,IEF,IP_F,nfmopcm
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      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
      DIMENSION mulfg(*),scffrg(*),fmoscf(*),fmoci(*),fmodft(*),fmocc(*)
     *         ,mpnfmo(*)
      data rhf/8HRHF     /,rnone/8HNONE    /,none/4HNONE    /,
     *     OPTMIZ /8HOPTIMIZE/,energy/8HENERGY  /,GRAD/8HGRADIENT/,
     *     OPTFMO/8HOPTFMO  /,rmc/8HMCSCF   /
c
      if(nder.gt.1) then
        if(maswrk) write(iw,8000)
        call abrt
      endif
      domcscf=.false.
      doci=.false.
      dodft=.false.
      docc=.false.
      domp=.false.
      legalscf=.true.
      legalmul=.true.
c 
      do ilay=1,nlayer
        if(fmoscf(ilay).eq.rmc) then
          domcscf=.true.
        else if(fmoscf(ilay).ne.rhf) then
          legalscf=.false.
        endif
        if(fmoci(ilay).ne.rnone) then
          doci=.true.
        endif
        if(fmodft(ilay).ne.rnone) then
          dodft=.true.
        endif
        if(fmocc(ilay).ne.rnone) then
          docc=.true.
        endif
        if(mpnfmo(ilay).ne.0) then
          domp=.true.
        endif
      enddo
      nmc=0
      do ifg=1,nfg
        if(scffrg(ifg).eq.rmc) nmc=nmc+1
        if(mulfg(ifg).ne.1.and.scffrg(ifg).ne.rmc) legalmul=.false. 
      enddo
      if(nmc.gt.1) legalscf=.false.
c
      if(nbody.gt.2.and.(RESPAP(1).ne.zero.or.rESPPC(1).ne.zero)) then
        if(maswrk) write(iw,9000)
        call abrt
      endif
      if(domcscf.and.nbody.gt.2) then
        if(maswrk) write(iw,9010)
        call abrt
      endif
      if(.not.legalscf) then
        if(maswrk) write(iw,9020)
        call abrt
      endif
      if(doci.and.nbody.ge.0) then
        if(maswrk) write(iw,9030)
        call abrt
      endif
      if(dodft.and.nbody.gt.3) then
        if(maswrk) write(iw,9040)
        call abrt
      endif
      if(docc.and.nbody.gt.3) then
        if(maswrk) write(iw,9050)
        call abrt
      endif
      if(domp.and.nbody.gt.2) then
        if(maswrk) write(iw,9060)
        call abrt
      endif
      if(runtyp.ne.energy.and.runtyp.ne.GRAD.and.runtyp.ne.OPTMIZ.and.
     *   runtyp.ne.OPTFMO) then
        if(maswrk) write(iw,9070)
        call abrt
      endif
      if(MPCTYP.ne.none) then
        if(maswrk) write(iw,9080)
        call abrt
      endif
      if(.not.legalmul) then
        if(maswrk) write(iw,9090)
        call abrt
      endif
      if(IECP.ne.0) then
        if(maswrk) write(iw,9100)
        call abrt
      endif
      if(nfmopcm.gt.nbody) then
c     if(nfmopcm.gt.0) then
        write(iw,9200)
        call abrt
      endif
      if(nfmopcm.ne.0.and.nder.gt.0) then
        write(iw,9210)
        call abrt
      endif
c
      return 
 8000 format(/1x,'FMO Hessian is not available.')
 9000 format(/1x,'RESPAP and RESPPC must be 0 for FMO3.',/)
 9010 format(/1x,'FMO3-MCSCF is not implemented.',/)
 9020 format(/1x,'Only SCFTYP=RHF or MCSCF is permitted.',/)
 9030 format(/1x,'FMO-CI is not implemented.',/)
 9040 format(/1x,'FMO4-DFT is not implemented.',/)
 9050 format(/1x,'FMO4-CC is not implemented.',/)
 9060 format(/1x,'FMO3-MP2 is not implemented.',/)
 9070 format(/1x,'Only RUNTYP=ENERGY, GRAD, OPTIMIZE or OPTFMO is ',
     *           'permitted.',/)
 9080 format(/1x,'Semiempirical runs not supported with FMO.',/)
 9090 format(/1x,'Only the MCSCF fragment may be non-singlet.',/) 
 9100 format(/1x,'No core potentials may be used with FMO.',/)
 9200 format(/1x,'IFMO in $PCM should not be exceed NBODY in $FMO.')
 9210 format(/1x,'FMO-PCM gradient is not available.')
      END
C*MODULE fmoprp   *DECK fmo0x
      SUBROUTINE fmo0x
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      Integer ddi_world,ddi_group
      Parameter(ddi_world=0,ddi_group=1)
      parameter(MXATM=2000,mxrt=100,maxelm=3,zero=0.0D+00,one=1.0D+00)
      logical isgddi,parout,GOPARR,DSKWRK,MASWRK,myjob,orbxch,enexch,
     *        docas,atonce,savdsk,densfmo0,dosap,dospc
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FMCOM / X(1)
      COMMON /GDDI/   ISCOPE,NGROUPS,MYGROUP,MEGLOB,npglob,nnglob,JBTYP,
     *                ISGDDI,PAROUT
      COMMON /GRDPAR/ ORIGIN(3),XVEC(3),YVEC(3),ZVEC(3),UX(3),UY(3),
     *                UZ(3),GRDSIZ,NGRID,IGUNIT,NXG,NYG,NZG,modgrd
      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 /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /XYZPRP/ XP,YP,ZP,
     *                DMX,DMY,DMZ,
     *                QXX,QYY,QZZ,QXY,QXZ,QYZ,
     *                QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ,
     *                OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                OXZZ,OYZZ,OZZZ,OXYZ,
     *                OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,convfg,cnvdmp,coroff,IXESP,mxitfg,nguess,
     *                NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c
      common /fmopnt/ lichfg,lmulfg,lidmrec,lfrgnam,llayfrg,lindat,
     *                lnCBS,lfmozan,lfmoc,lizbas,liaglob,libdgh,liabdfg,
     *                ljabdfg,lnCAO,lidxCAO,liaprjo,ljaprjo,lCoreAO,
     *                lOccCor,lshiftb,liodfmo,lfmoda,lfmodb,lfmoespa,
     *                lfmoespb,llocfmo,lscffrg,lrij,lpopmul,lpopmat,
     *                lialoc,lindbd,liatfrg,lindfrg,lindgfrg,lnatfrg,
     *                lnat0frg,lianfrg,lzanfrg,lcfrg,llibish,llibnsh,
     *                llibng,lindatg,lfmobuf(2),lfmode,lnumfrg,lloctat,
     *                liaoglob,lloadm,lfmoge,ldgrid,liodcfmo
      common /fmorun/ espscf,IDAFMO,icurlay,icurfg,jcurfg,kcurfg,ncursh,
     *                ngau,icurpop,ifmostp,moncor,nat1e,needr,modrst,
     *                norbproj
      dimension nelm(0:maxelm)
      data nelm/0,3,9,19/
c
c     Generate the proper 0-body FMO wavefunction.
c     0-body means free monomers. If fractioned bonds are present, this 
c     means the monomer state while capped with methyls.
c     In other words, for molecule F1-F2-F3 the 0-body state is got as:
c     F1'+F2'+F3', where
c     F1' is in F1'-Me 
c     F2' is in Me-F2'-Me 
c     F3' is in    F3'-Me 
c     Me means methyl, so the fractioned bonds are C-C. 
c     The code is general for other bonds, just more library data should be
c     entered. For example, if C-N bond is fractioned, the cap is NH2.
c
c     The calculation proceeds by calling the FMO driver for each monomer.
c     Capping has to be done in FMOX where all information is stored.
c
      if(maswrk) write(iw,9000) 
c
c     Call fmox to read in most options.
c
      call fmox(0)
c
      nfmoelm=0
c
c     Allocate memory for non-sparse grid. It should be processed here to
c     minimise the overhead. 
c
      nxyzg=0
      if(iand(modprp,4).ne.0) nxyzg=nzg*nyg*nxg
      nfmoelm=nelm(iemom)
      nprocg=nproc
      ngrmax=nprocg
      orbxch=mod(modorb,2).ne.0
      enexch=mod(modorb/2,2).ne.0
      atonce=iand(modpar,4).ne.0
      densfmo0=maxl1c.ne.0
      dosap=respap(1).ne.zero.or.respap(2).ne.zero
      dospc=resppc(1).ne.zero.or.resppc(2).ne.zero
c
cnb   This is a despicable way to set DOCAS.
c
      docas=orbxch.and.enexch
c
      maxl2c=(maxl1c*maxl1c+maxl1c)/2
      maxl3c=maxl1c*maxl1c
      n30ex=5+nlayer
      LDAR30c=maxl2c
      if(orbxch) LDAR30c=maxl3c
      if(docas)  LDAR30c=maxl2c+maxl3c
      if(enexch) LDAR30c=LDAR30c+maxl1c
      NDAR30c=1+nfg*2*2
      NDAR30ci=(NDAR30c+n30ex-1)/nwdvar+1
      LDAR30c=max(LDAR30c,NDAR30ci)
      IDAcFMO=30
c
      maxm2d=LDAR30c
      if(dosap) maxm2d=maxm2d+maxl1c
      if(dospc) maxm2d=maxm2d+maxl1c
c
      CALL VALFM(LOADFM)
      memon=LOADFM+1
      mfmoelm=memon+nfg*4
      mfrgnam=mfmoelm+nfmoelm*(nfg+1)
      mlayfrg=mfrgnam+nfg
      ldgrid=mlayfrg+(nfg-1)/nwdvar+1
      mjobgrp=ldgrid+nxyzg
      mnumfrg=mjobgrp+(nfg+1-1)/nwdvar+1
      mnatfrg=mnumfrg+(nfg-1)/nwdvar+1
      mscffrg=mnatfrg+(nfg-1)/nwdvar+1
      liodcfmo=mscffrg+nfg
      mmaxl30=liodcfmo+NDAR30ci
      mmastid=mmaxl30+(nfg-1)/nwdvar+1
      mfrgscf=mmastid+(nprocg-1)/nwdvar+1
      last=mfrgscf+nfg
      mfmoda=last
      mifgbuf=last
      if(atonce) then
        mnfgbuf=mifgbuf+(2*nfg*ngrmax-1)/nwdvar+1
        last=mnfgbuf+(2*ngrmax-1)/nwdvar+1
      else
        mnfgbuf=last
        last=mfmoda+maxm2d
c       lfmoda is surreptitiosly used in dmexch (if not atonce).
c       Note that lfmoda is reset in FMOX, so it has to be set after
c       calling it.
      endif
c
      if(nxyzg.ne.0) then
c       It is possible to save some memory in FMOX if the arrays are allocated 
c       below. 
        mfmozan=last
        mfmoc=mfmozan+natfmo
        mizbas=mfmoc+3*(natfmo+1)
        last=mizbas+(natfmo-1)/nwdvar+1
      endif
      NEED = LAST- LOADFM -1
      CALL GETFM(NEED)
      if(maswrk) write(iw,9010) NEED
      call vclr(x(memon),1,nfg*4)
      call vclr(x(mfmoelm),1,nfmoelm*(nfg+1))
      call viclr(x(mlayfrg),1,nfg)
      call vclr(x(mfrgnam),1,nfg)
      if(nxyzg.ne.0) call vclr(x(ldgrid),1,nxyzg)
      call viclr(x(mjobgrp),1,nfg)
      call viclr(x(mnumfrg),1,nfg)
      call viclr(x(mnatfrg),1,nfg)
      call vclr(x(mscffrg),1,nfg)
      call viclr(x(mmaxl30),1,nfg)
c
      savdsk=DSKWRK
      if(densfmo0) then
        call viclr(x(liodcfmo),1,NDAR30c+n30ex)
        DSKWRK=.true.
        call RAOPEN2(IDAcFMO,x(liodcfmo),0,NDAR30c,LDAR30c,0,NPRINT)
c       It is vital to reserve the first record on IDAcFMO for indexing.
        call rawrites(IDAcFMO,x(liodcfmo),x(liodcfmo),NDAR30ci,NDAR30ci,
     *                1,0)
        DSKWRK=savdsk
      endif
c
      if(isgddi) then
        call gddi_scope(ddi_group)
        call gddi_mastid(x(mmastid))
        call GDDICOUNT(-1,lgroup,myjob)
      else
        call ixstor(x(mmastid),1,master)
c       mastid(0)=master
      endif
c
      do 200 ifg=1,nfg
        if(isgddi) then
          call GDDICOUNT(0,lgroup,myjob)
          if(.not.myjob) goto 200
        endif
        call fmox(ifg)
        call ixstor(x(mjobgrp),ifg,mygroup)
        if(maswrk) then
          x(memon+ifg-1  )=ETOT 
          x(memon+ifg-1+nfg)=ESCF
          if(nfmoelm.ne.0) 
     *      call dcopy(nfmoelm,DMX,1,x(mfmoelm+(ifg-1)*nfmoelm),1)
          call ixstor(x(mlayfrg),ifg,icurlay)
          call ixstor(x(mnumfrg),ifg,num)
          call ixstor(x(mnatfrg),ifg,nat)
          call ixstor(x(mmaxl30),ifg,nat1e)
          x(mfrgnam+ifg-1)=espscf
          x(mfrgscf+ifg-1)=scftyp
        endif
  200 continue
c
      if(isgddi) then
        call GDDICOUNT( 1,lgroup,myjob)
        call gddi_scope(ddi_world)
c
c       The results are harvested from all GDDI masters (in the world scope).
c
        if(nfmoelm.gt.0) call ddi_gsumf(2419,x(mfmoelm),nfmoelm*nfg)
        call ddi_gsumi(2419,x(mlayfrg),nfg)
        call ddi_gsumi(2419,x(mnumfrg),nfg)
        call ddi_gsumi(2419,x(mnatfrg),nfg)
        call ddi_gsumi(2419,x(mmaxl30),nfg)
c
c       The operation below is potentially dangerous, since we are collecting 
c       Hollerith characters in real precision. But the probability of 
c       hitting something like NAN is nearly zero (yet the wise say that a cow
c       drowned where it was knee-deep on the average).
c
        call ddi_gsumf(2419,x(mfrgnam),nfg)
        call ddi_gsumf(2419,x(mfrgscf),nfg)
c
c     Exchange the density file
c
        if(densfmo0) then
          call gddi_scope(ddi_group)
c         Some thought is needed to see if this will work in the multilayer
c         case. nconv is set to -1 to prevent summing of mglocon, there being
c         no space allocated for it).
c         dmexch needs to be called in the group scope after indexing
c         arrays have been summed in the world scope above. 
          ilay=1
          nfge=nfg*2
          nconv=-1
          irec0=1
          lfmoda=mfmoda
          call dmexch(ilay,nfge,x(mjobgrp),x(mlayfrg),x(mnumfrg),
     *                x(mnatfrg),x(mscffrg),IDAcFMO,x(liodcfmo),
     *                x(mmaxl30),x(memon),dum,nconv,irec0,.false.,
     *                .false.,orbxch,enexch,atonce,x(mmastid),dum,dum,
     *                x(mifgbuf),x(mnfgbuf))
c         dum's are mglocon and orb+at populations that need not be exchanged. 
          call gddi_scope(ddi_world)
        else
c         In this case memon has not been summed in dmexch.
          call ddi_gsumf(2419,x(memon),nfg*4)
        endif
      endif
c
c     Output the results (cloned from fmoprop).
c
      if(maswrk) then  
        write(iw,9114)
        nenm=1
        if(rcorsd.ne.0) nenm=2
        if(nenm.gt.1) write(iw,9115)
        if(nenm.le.1.and.iemom.ge.1) write(iw,9116)
        emontot=zero
        emontotu=zero
        mfmoelmt=mfmoelm+nfg*nfmoelm-1
        do ifg=1,nfg
          nelm1=nelm(1)
          if(iemom.lt.1) nelm1=0
          eifg= x(memon+ifg-1  )
          eifgu=x(memon+ifg-1+nfg)
          ilay=ixftch(x(mlayfrg),ifg)
          frgnami=x(mfrgnam+ifg-1)
          mfmoelmi=mfmoelm+(ifg-1)*nfmoelm-1
          call daxpy(nfmoelm,one,x(mfmoelmi+1),1,x(mfmoelmt+1),1)
          if(nenm.eq.2) then
            write(iw,9118) ifg,frgnami,ilay,eifg,eifgu,
     *                     (x(mfmoelmi+ielm),ielm=nelm(0)+1,nelm1)
          else
            write(iw,9120) ifg,frgnami,ilay,eifg,
     *                    (x(mfmoelmi+ielm),ielm=nelm(0)+1,nelm1) 
          endif
          if(iemom.ge.2)
     *      write(iw,8122) (x(mfmoelmi+ielm),ielm=nelm(1)+1,nelm(2))
          if(iemom.ge.3)
     *      write(iw,8123) (x(mfmoelmi+ielm),ielm=nelm(2)+1,nelm(3))
          emontot=emontot+eifg
          emontotu=emontotu+eifgu
        enddo
        if(nenm.ge.2) then
          write(iw,9200) '  corr',0,emontot
          write(iw,9200) 'uncorr',0,emontotu
          write(iw,9200) 'delta ',0,emontot-emontotu
        else
          write(iw,9200) 'uncorr',0,emontot
        endif
        if(iemom.ge.1)
     *    write(iw,9230) 0,(x(mfmoelmt+ielm),ielm=nelm(0)+1,nelm(1)),
     *                   DNRM2(3,x(mfmoelmt+1),1)
        if(iemom.ge.2) 
     *    write(iw,9231) (x(mfmoelmt+ielm),ielm=nelm(1)+1,nelm(2))
        if(iemom.ge.3) 
     *    write(iw,9232) (x(mfmoelmt+ielm),ielm=nelm(2)+1,nelm(3))
c
        etot=emontot
        write(iw,9900) 
      endif
c
c     Write up a nice restart data group
c
      if(nlayer.eq.1) then
        ilay=1
        call eminout(1,ilay,0,x(mlayfrg),x(memon),.false.) 
      endif
c
      if(densfmo0) then
c       l1fmo should really be computed in FMO0 but until then just reset to 0.
        idum=0
        l1fmo=0
        l1fmoc=0
        DSKWRK=.true.
        call bconvrec(dum,idum,idum,NDAR30c,x(liodcfmo),l1fmo,l1fmoc,
     *                IDAcFMO,idum,.false.)
        call rawrites(IDAcFMO,x(liodcfmo),x(liodcfmo),NDAR30ci,NDAR30ci,
     *                1,0)
        CALL RACLOS(IDAcFMO,'KEEP')
        DSKWRK=savdsk
      endif
c
c     Finally finish off by writing the grid data (if any).
c
      if(nxyzg.ne.0) then
c       grid density has to be summed for both GDDI and regular DDI.
        call ddi_gsumf(2416,x(ldgrid),nxyzg)
c       Since FMOX destroys coordinates+charges, they have to be rearead.
        call fmoxyz(' $FMOXYZ',0,natfmo,x(mfmozan),x(mfmoc),x(mizbas))
        if(maswrk)
     *    call makecube(0,0,natfmo,x(mfmozan),x(mfmoc),x(ldgrid),
     *                  NXG,NYG,NZG,ORIGIN,UX,UY,UZ)
      endif
      CALL RETFM(need)
      call timit(1)
      return
 9000 format(/1x,76(1H=),
     *       /31x,'Zero-body FMO driver.',
     *       /17x,'The purpose is to obtain free monomer properties.'
     *       /1x,'For covalently bonded fragments free monomers ',
     *           'are obtained by adding capping',
     *       /10x,'fragments (such as methyl) at the fractioned bonds.',
     *       /1x,76(1H=)/)
 9010 format(/1x,'Using',I10,' words of memory.',/)
 9114 format(/1x,'Zero-body FMO properties.',/1x,25(1H=)/)
 9115 format(26x,'Ecorr',11x,'Euncorr',7x,'DX',7x,'DY',7x,'DZ')
 9116 format(28x,'E',11x,'DX',7x,'DY',7x,'DZ')
 9118 format(1x,I4,'(',A8,',L',I1,')',2F17.9,3F9.5)
 9120 format(1x,I4,'(',A8,',L',I1,')', F17.9,3F9.5)
 8122 format(8x,'Q=',6F10.5)
 8123 format(8x,'O=',5F10.5,/8x,2x,5F10.5)
 9200 format(/1x,'Total energy of the molecule: E',A6,'(',I1,')=',F19.9)
 9230 format(1x,'Dipole moment D(xyz),DA(',I1,')=',4F13.7)
 9231 format(1x,'Quadrupole moment Q(XX,YY,ZZ)=',3F12.7/,
     *       1x,'                   (XY,XZ,YZ)=',3F12.7)
 9232 format(1x,'Octopole moment O(XXX,XXY,XXZ,XYY)=',4F11.6/,
     *       1x,'                 (YYY,YYZ,XZZ,YZZ)=',4F11.6/,
     *       1x,'                         (ZZZ,XYZ)=',2F11.6)
 9900 format(/)
      end
c
C*MODULE fmoprp  *DECK fmototd 
      SUBROUTINE fmototd(ida,l2,da,totd,iaoglob)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION da(*),totd(*),iaoglob(*)
      logical GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
c
c     Accumulate the electron density for the total system in AO basis.
c     parstat: GroupFull (statically).
c
      if(ida.eq.0) return
c
      nsize=(l2-1)/nproc+1
      myst=me*nsize+1
      myfi=min(myst+nsize-1,l2)
      call tribrk(myst,-1,i,j)
      i=i-1
      fda=ida
c
      do ij=myst,myfi
        ig=iaoglob(i) 
        jg=iaoglob(j)
c       write(6,*) 'wwwij',ij,i,j,ig,jg
        if(jg.lt.ig) then
          ijg=(ig*ig-ig)/2+jg
        else
          ijg=(jg*jg-jg)/2+ig
        endif
        totd(ijg)=totd(ijg)+fda*da(ij)
        j=j+1
        if(j.gt.i) then
          j=1
          i=i+1
        endif
      enddo
c
      return
      end
c
C*MODULE fmoprp  *DECK addecc
      subroutine addecc(ibody,ida1,ida2,ida3,extracc,nextracc)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL CNVR,QDR,QDCR
      COMMON /CCENGY/ ENRG,EREF,EMP2,eccn,ETOT(6),ECORR(6),
     *                DIAGS(3),AMPMX(5,2),IAMPMX(5,4,2),XO1,XO2,
     *                DIFMAX,DIFFENG,ITER,CNVR
      COMMON/QUADRUPLE/QDR,QDCR
      COMMON/T4TOT/ETT2,ETS2,EQ1A,EQ1B,EQ2A,EQ2B,EQ3A,EQ3B,EQ4A,
     *EQ4B,O4A,O4B,ECORQ(8),ETOTQ(8)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      dimension extracc(nextracc,*),fda(3)
c
c     Accumulate extra CC properties.
c
      fda(1)=ida1
      fda(2)=ida1+ida2
      fda(3)=ida1+ida2+ida3
      do ii=ibody,nbody
        fdai=fda(ii)
        if(fdai.ne.0.0D+00) then
          extracc(1,ii)=extracc(1,ii)+fdai*emp2
          extracc(2,ii)=extracc(2,ii)+fdai*eccn
          call daxpy(6,fdai,ecorr,1,extracc(3,ii),1)
          if(QDR) call daxpy(7,fdai,ECORQ(2),1,extracc(9,ii),1)
        endif
      enddo
c
      return
      END
c
C*MODULE fmoprp  *DECK pcmnup
      SUBROUTINE pcmnup(ida,vnuc)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      logical GOPARR,DSKWRK,MASWRK
      PARAMETER (MXATM=2000,zero=0.0D+00)
      dimension vnuc(*)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FE,UNZOR,PORT
      COMMON /PCMDIM/ mxsp,mxts,mempcm1,mempcm2,NTS
      COMMON /PCMPNT/ LXYZRE,LSSFE,LLIST,LALPCM,LRINPCM,LINAPCM,LINFPCM,
     *                LINIPCM,LMEPCM,LAXYZCT,LXYZCT2,LNVERT,LQSN,LQSE,
     *                LQFS,LQIND,LISPHE,LVAD,LQOR,LVECMUL,LAIPRJ,LFIPRJ
     *                ,LPEL
C
C     Compute nuclear potential due to induced charges in the FMO-PCM.
C
c     if(ida.eq.0) return
      da=ida
c     write(6,*) 'wwwnup',ida
C
C     THE NUCLEAR POTENTIAL GOES TO VNUC
c     PC will contain nuclear potential contribution to the n-mer energy.
c     VNUC accumulates potential for the whole system, thus it has the DA
c     prefactor, whereas PC is computed for the n-mer only, no prefactor.
c     VNUC is not globally summed here, but PC is.
c
cnb   The present scheme is extravagant. The right way is to compute PCs and
c     VNUC for monomers only, since dimer values are additive of monomers! 
C
      PC=zero
      DO 100 ITS=1,NTS
        if(goparr.and.MOD(ITS,NPROC).NE.me) goto 100
        VITS=zero
        xctsi=x(LAXYZCT+ITS-1)
        yctsi=x(LAXYZCT+MXTS+ITS-1)
        zctsi=x(LAXYZCT+MXTS*2+ITS-1)
        DO IAT=1,NAT
          VITS=VITS+ZAN(IAT)/SQRT((XCTSI-C(1,IAT))**2+
     *                       (YCTSI-C(2,IAT))**2+(ZCTSI-C(3,IAT))**2)
        ENDDO
        if(ida.ne.0) VNUC(ITS)=VNUC(ITS)+da*VITS
        PC=PC+VITS*x(lQSE+ITS-1)
  100 continue
      if(goparr) call ddi_gsumf(2416,PC,1)
      RETURN
      END
C*MODULE fmoprp  *DECK calcasc
      SUBROUTINE calcasc(itermscf,ntspar,ichfmo,nefmo,lfmasc,needasc,
     *                   lesolv,nesolv,lpel,lvnuc)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c
      logical GOPARR,DSKWRK,MASWRK,PACK2E
c 
      PARAMETER (ONE=1.0D+00,TOANGS=0.52917724924D+00,ANTOAU=one/TOANGS)
c
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      common /INT2IC/ NINTIC,iNINTIC,nxxic,lbufpic,lixic,labsix,nintix
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /PCKLAB/ LABSIZ
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FE,UNZOR,PORT
      COMMON /PCMDIM/ mxsp,mxts,mempcm1,mempcm2,NTS
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IRPPCM,IEF,IP_F,nfmopcm
      COMMON /PCMPLY/ STOT,VOL,CCX,CCY,CCZ,RDIF,NESF,NESFP,NESFF,I_NESF,
     *                L_AST
      COMMON /PCMPNT/ LXYZRE,LSSFE,LLIST,LALPCM,LRINPCM,LINAPCM,LINFPCM,
     *                LINIPCM,LMEPCM,LAXYZCT,LXYZCT2,LNVERT,LQSN,LQSE,
     *                LQFS,LQIND,LISPHE,LVAD,LQOR,LVECMUL,LAIPRJ,LFIPRJ
     *                ,LPELDUM
c
c     Compute apparent surface charges (ASC). 
c     itermscf=0 estimate memory need 
c     itermscf=1 do the work 
c     No memory is allocated inside.
c
C     -- MEMORY FOR SALVATION
C
      MXDII1 = MXDIIS + 1
c
      LCHG   = lfmasc + 1
      LQPOT  = LCHG   + NTS
      LQ0    = LQPOT  + NTS
      LQ1    = LQ0    + NTS
      LQ2    = LQ1    + NTS
      LQ3    = LQ2    + NTS
      LD0    = LQ3    + NTS
      LQA    = LD0    + NTS
      LDIMAT = LQA    + NTS
      LTMP0  = LDIMAT + MXDII1*MXDII1
      LTMP1  = LTMP0  + NTS*3
      LTMP2  = LTMP1  + MXDII1
      LIPVT  = LTMP2  + MXDII1*MXDII1
      LRMUL  = LIPVT  + MXDII1
      lesolv = LRMUL  + MXSP*10
      lpel   = lesolv + nesolv
      lvnuc  = lpel   + NTS
      last   = lvnuc  + NTS
c
c     Allocate the largest PCM array.
c     Try borrowing integral memory, if available.
c     Note that OPTFMOX also uses it. 
c     We can use the integral buffer here since QREP is used in ASCIT
c     after SCF is all done and the integrals are no longer needed. 
c     Try allocating all PCM in integrals?
c
      nqrep=NTSPAR*MXDII1*2
      iuse2e=0
      if(NINTIC.gt.0) then
        LABSIZ2=2/LABSIZ
        n2ebuf=NINTMX+NINTIC+(NINTMX+NINTIC-1)/LABSIZ2+1
c       write(6,*) 'Buffer',n2ebuf,nhblock
        if(n2ebuf.ge.nqrep) then
          LQREP=LBUFPIC
          iuse2e=1
        endif
      endif
      if(iuse2e.eq.0) then
        LQREP = last 
        last  = LQREP  + nqrep
      endif
c
      NEEDasc = LAST -lfmasc -1
      if(itermscf.eq.0) return
c
      if(maswrk) then
        write(iw,9000) NEEDasc
        if(iuse2e.ne.0) write(iw,9010) nqrep
      endif
c
c     Add electron and nuclear potentials and normalise.
c
c     write(6,*) 'wwwpcm-elpot'
c     call prsqa(x(lPEL),NTS,1,1)
c     call abrt
c
c     For charge compensation, solve PCM equations first for nuclear
c     potential/charge only, then full equations.
      if(icomp.ne.0.and.itermscf.eq.1) then
c       Nuclear charges are got only once in case of charge compensation.
c       Try to give more iterations to ensure convergence.
        call dcopy(NTS,x(lVNUC),1,X(LQPOT),1)
        call dscal(NTS,-one,X(LQPOT),1)
        MXIsav1=MXITR1
        MXIsav2=MXITR2
        MXITR1=max(mxitr1,2000)
        MXITR2=max(mxitr2,2000)
        CALL ASCIT(itermscf,X(LQPOT),X(LQ0),X(LQ1),X(LQ2),X(LQ3),X(LD0),
     *             X(LQA),X(LDIMAT),X(LQREP),X(LTMP0),X(LTMP1),X(LTMP2),
     *             X(LIPVT),X(LRMUL),MXDII1,ntspar,x(LAXYZCT),
     *             x(LAXYZCT+MXTS),x(LAXYZCT+MXTS*2),x(LAXYZCT+MXTS*3),
     *             x(LQSE),x(LXYZRE),x(LXYZRE+MXSP),x(LXYZRE+MXSP*2),
     *             x(LXYZRE+MXSP*3),x(LISPHE),x(LLIST))
        MXITR1=MXIsav1
        MXITR2=MXIsav2
        IF(IEF.EQ.10) call dscal(NTS,(EPS-ONE)/EPS,x(LQSE),1)
        call dcopy(NTS,x(LQSE),1,x(LQSN),1)
        call timit(1)
        if(maswrk) write(iw,9050)
      endif
      if(icomp.ne.0) then 
c       Store just the electron potential in QPOT for the next call to ASCIT.
        call dcopy(NTS,x(lPEL),1,X(LQPOT),1)
        call dscal(NTS,-one,X(LQPOT),1)
      else 
c       Store the total potential in QPOT for the next call to ASCIT.
        call vadd(x(lpel),1,x(lVNUC),1,X(LQPOT),1,NTS)
        call dscal(NTS,-one,X(LQPOT),1)
      endif
c
c     write(6,*) 'wwwpcm-nucpot unnorm'
c     call prsqa(x(lVNUC),NTS,1,1)
c     write(6,7777) (X(LQPOT+i),i=0,NTS-1)
c7777 format(1x,5F15.8)
c     call prsqa(X(LQPOT),NTS,1,1)
C
      CALL ASCIT(itermscf,X(LQPOT),X(LQ0),X(LQ1),X(LQ2),X(LQ3),X(LD0),
     *           X(LQA),X(LDIMAT),X(LQREP),X(LTMP0),X(LTMP1),X(LTMP2),
     *           X(LIPVT),X(LRMUL),MXDII1,ntspar,x(LAXYZCT),
     *           x(LAXYZCT+MXTS),x(LAXYZCT+MXTS*2),x(LAXYZCT+MXTS*3),
     *           x(LQSE),x(LXYZRE),x(LXYZRE+MXSP),x(LXYZRE+MXSP*2),
     *           x(LXYZRE+MXSP*3),x(LISPHE),x(LLIST))
C
C  SCALE C-PCM CHARGES BY (EPS-1)/EPS
C
      IF(IEF.EQ.10) call dscal(NTS,(EPS-ONE)/EPS,x(LQSE),1) 
c
c     write(6,*) 'wwwpcm-asc'
c     call prsqa(QSE,NTS,1,1)
c     ENPCM=-ddot(NTS,X(LQPOT),1,QSE,1)
c     px=ddot(NTS,x(lPEL),1,QSE,1)
C     PB = INTERACTION ELECTRONS-NUCLEAR INDUCED CHARGES
C     PC = INTERACTION NUCLEI-ELECTRONIC INDUCED CHARGES
C     PX = INTERACTION ELECTRONS-ELECTRONIC INDUCED CHARGES
c     PB=zero
c     PC=-px
c     PX=ENPCM - PX
c
c     FMO-PCM differs from the ab initio PCM in putting ASC-electron integrals.
c     In ab initio they go to the Fock matrix directly.
c     In FMO-PCM they are added to one-electron integrals (and thus to
c     the Fock matrix).
c     The Fock matrices are the same but the 1e inetgrals differ. Thefore, the
c     energy correction term differs (subtracting 1/2*Tr(D*Vel)=1/2*(Pel*Qse)).
c     P.S. The reason for the different behaviour is that in FMO-PCM the ASCs
c     are constant in SCF iterations.
C
C     --- NORMALIZATION ---
C
C  4) NORMALIZZAZIONE DELLE CARICHE VIRTUALI ELEMENTARI ELETTRONICHE.
C     PER ICOMP=1 LA CORREZIONE E' DISTRIBUITA IN PROPORZIONE
C     ALL'AREA DELLA TESSERA.
C     PER ICOMP=2 LA CORREZIONE E' EFFETTUATA CON FATTORI COSTANTI.
C
      if(icomp.ne.0) then
        QN=-(nefmo+ichfmo)*(EPS-ONE)/EPS
        QE=nefmo*(EPS-ONE)/EPS
        QTE=ddot(NTS,x(LQSE),1,one,0)
        QTN=ddot(NTS,x(LQSN),1,one,0)
c       QSE=Qe, QSN=Qn 
c       write(6,*) 'wwwCC',QE,QTE,QN,QTN
        IF(ICOMP.EQ.1) THEN
c       charge renormalisation is:
c       Qt=Qt+Fe*S+Fn*S=Qe+Qn+(Fn+Fe)*S
          SUPTOT=STOT*ANTOAU*ANTOAU
          FE=(QE-QTE)/SUPTOT
          FN=(QN-QTN)/SUPTOT
          call daxpy(NTS,ONE,x(LQSN),1,x(LQSE),1)
c         x(LAXYZCT+MXTS*3) is AS
          call daxpy(NTS,FN+FE,x(LAXYZCT+MXTS*3),1,x(LQSE),1)
        ELSE IF(ICOMP.EQ.2) THEN
c       charge renormalisation is:
c       Qt=Qe*Fe+Qn*Fn
          FE=QE/QTE
          FN=QN/QTN
          call dscal(NTS,FE,x(LQSE),1)
          call daxpy(NTS,FN,x(LQSN),1,x(LQSE),1)
        END IF
        QETN=ddot(NTS,x(LQSE),1,one,0)
        if(maswrk) WRITE(IW,9200) QTE+QTN,QETN,FN,FE
c       Note that FE and FN have very different meanings in ICOMP=1 and 2.  
      else
c       Without charge compensation, QSE contains Qe+Qn already. 
        QET=ddot(NTS,x(LQSE),1,one,0)
        TCH=-ichfmo*(EPS-ONE)/EPS
        if(maswrk) WRITE(IW,9100) QET,TCH
      END IF
c
c     We must be in world scope now, so broadcast charges to ensure
c     no SCF deadlocks. 
c
c     if(goparr) CALL DDI_BCAST(2422,'F',QSE,NTS,0)
c
      call timit(1)
c
      RETURN
 9000 format(/1x,'Calculating apparent surface charges...',
     *       /5x,'Using',I12,' words of memory.')
 9010 format(1x,I12,' words are borrowed from the integral buffer.')
 9050 format(1x,'Calculated nuclear charges, now doing electronic...')
 9100 format(1x,'PCM ASC=',F18.10,' theoretical=  ',F18.10,/)
 9200 format(1x,'PCM ASC=',F18.10,' normalised to=',F18.10,/
     *      1x,'Nuclear factor=',F15.10,', electronic factor=',F15.10,/)
      END
