      subroutine scfcycle(iteration,nreal,ninteger,nmaxatoms,nmaxorbs,
     &                    maxaux,nmaxpoints,idiis,nfunctional,ngridtype,
     &                    natoms,nptcharges,nbsystems,level,nalpha,
     &                    nbeta,ncontractions,nconts,ncontp,ncontd,ncds,
     &                    ncdspd,nxcs,nxcspd,ixcfit,natomtype,nptsatom,
     &                    nshels,nshelp,nsheld,mtloca,ilfunc,icfunc,
     &                    ngaussians,nfuncatom,icdlfunc,icdcfunc,
     &                    ixclfunc,ixccfunc,isys,icdlsave,iwkvec,energy,
     &                    eshift,dmixing,diiserror,dmaxpc,cutoff,coord,
     &                    coordpc,charge,ptcharge,core,focka,fockb,
     &                    focke,dmta,dmtb,coeffa,coeffb,overlap,alpha,
     &                    coeffs,alphacd,coefscd,coefpcd,coefdcd,
     &                    alphaxc,coefsxc,coefpxc,coefdxc,cdfitc,xcfita,
     &                    xcfitb,xcfite,tvector,tveca,tvecb,tvece,
     &                    storecd,cdstore,wkvec)

c  this subroutine written by alain st-amant of the department of
c  pharmaceutical chemistry, university of california, san francisco.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      include "mpif.h"

      dimension natomtype(*),nptsatom(*),nshels(*),nshelp(*),nsheld(*),
     &          mtloca(*),ilfunc(*),icfunc(*),ngaussians(*),
     &          nfuncatom(*),icdlfunc(*),icdcfunc(*),ixclfunc(*),
     &          ixccfunc(*),isys(nmaxatoms,*),icdlsave(*),iwkvec(*)

      dimension cutoff(*),coord(3,*),coordpc(3,*),charge(*),ptcharge(*),
     &          core(*),focka(*),fockb(*),focke(*),dmta(*),dmtb(*),
     &          coeffa(nmaxorbs,*),coeffb(nmaxorbs,*),
     &          overlap(nmaxorbs,*),alpha(*),coeffs(*),alphacd(*),
     &          coefscd(*),coefpcd(*),coefdcd(*),alphaxc(*),coefsxc(*),
     &          coefpxc(*),coefdxc(*),cdfitc(*),xcfita(*),xcfitb(*),
     &          xcfite(*),tvector(*),tveca(*),tvecb(*),tvece(*),
     &          storecd(*),cdstore(*),wkvec(*)

      data zero,one/0.0d0,1.0d0/

      call mpi_comm_rank(mpi_comm_world,myid,ierr)
      call mpi_comm_size(mpi_comm_world,numprocs,ierr)

      ndim=ncontractions*(ncontractions+1)/2

      ncdfunc=ncds+10*ncdspd
      nxcfunc=nxcs+10*nxcspd

      if(iteration.eq.1) then

               call startscf(nreal,ninteger,nmaxatoms,nmaxorbs,ncdfunc,
     &               nmaxpoints,natoms,nptcharges,nconts,ncontp,ncontd,
     &               ncds,ncdspd,nxcs,nxcspd,level,ngridtype,nbsystems,
     &               isys,icdlsave,iwkvec,nshels,nshelp,nsheld,ilfunc,
     &               icfunc,ngaussians,mtloca,icdlfunc,icdcfunc,
     &               ixclfunc,ixccfunc,natomtype,nptsatom,energy,dmaxpc,
     &               diiserror,cutoff,wkvec,core,focka,fockb,focke,
     &               overlap,coeffa,coeffb,coord,coordpc,charge,
     &               ptcharge,alpha,coeffs,alphacd,coefscd,coefpcd,
     &               coefdcd,alphaxc,coefsxc,coefpxc,coefdxc,cdfitc,
     &               xcfita,xcfitb,xcfite,dmta,dmtb,tveca,tvecb,tvece)

               call mpi_bcast(nptsatom,natoms,
     &                        mpi_integer,0,mpi_comm_world,ierr)

               if(nxcs.eq.0) nxcspd=0
               if(nxcs.eq.0) nxcfunc=0

               nnn=max(nxcfunc,ncontractions)

               ja=1
               jb=ja+3*nnn
               jc=jb+3*nnn
               jd=jc+nnn*nnn

               call gridcall(nreal-jd,nnn,1,iteration,nfunctional,
     &                       natoms,nmaxatoms,nbsystems,nalpha,nbeta,
     &                       nconts,ncontp,ncontd,nxcs,nxcs+nxcspd,
     &                       nmaxpoints,natomtype,nptsatom,nshels,
     &                       nshelp,nsheld,ilfunc,icfunc,ngaussians,
     &                       ixclfunc,ixccfunc,isys,iwkvec,
     &                       iwkvec(nnn+1),dmixing,cutoff(2),coord,
     &                       alpha,coeffs,alphaxc,coefsxc,coefpxc,
     &                       coefdxc,xcfita,xcfitb,xcfite,dmta,dmtb,
     &                       tveca,tvecb,tvece,wkvec(ja),wkvec(jb),
     &                       wkvec(jc),focka,fockb,focke,wkvec(jd))

                         endif

      call twoemm(nreal-3*ndim,ninteger,nconts,ncontp,ncontd,ncds,
     &            ncdspd,nxcs,nxcspd,nshels,nshelp,nsheld,mtloca,ilfunc,
     &            ngaussians,icfunc,icdlfunc,icdcfunc,ixclfunc,ixccfunc,
     &            iwkvec,cutoff(3),coord,alpha,coeffs,alphacd,coefscd,
     &            coefpcd,coefdcd,alphaxc,coefsxc,coefpxc,coefdxc,
     &            cdfitc,xcfita,xcfitb,xcfite,wkvec(1),focka,fockb,
     &            focke,wkvec(ndim+1),wkvec(2*ndim+1))

      if(ixcfit.eq.0) then
                            nxcs=0
                            nxcspd=0
                            nxcfunc=0
                      endif

      do 1001 j=1,ndim
      focka(j)=focka(j)+core(j)
      fockb(j)=fockb(j)+core(j)
 1001 focke(j)=focke(j)+core(j)

      ja=1
      jb=ja+ncontractions*ncontractions

      if((eshift.ne.zero).
     &    and.(myid.eq.0)) call carbo(nmaxorbs,ncontractions,
     &                                ncontractions,nalpha,nbeta,eshift,
     &                                focka,fockb,coeffa,coeffb,overlap,
     &                                wkvec(ja),wkvec(jb))

      ka=1
      kb=ka+ncontractions*ncontractions
      kc=kb+ncontractions*ncontractions
      kd=kc+ncontractions*ncontractions
      ke=kd+ncontractions*ncontractions
      kf=ke+ncontractions*ncontractions
      kg=kf+ncontractions*ncontractions
      kh=kg+ncontractions

      if(idiis.eq.1) then
                           if(kh.gt.nreal) stop 'calling diis'

                           call ddiis(81,iteration,ncontractions,
     &                                dmixing,diiserror,focka,fockb,
     &                                dmta,dmtb,wkvec(ka),wkvec(kb),
     &                                wkvec(kc),wkvec(kd),wkvec(ke),
     &                                wkvec(kf),wkvec(kg))
                     endif

      if(iteration.eq.1) tol=one
      if(iteration.gt.1) tol=dmaxpc**2

      dmaxpc=zero

      ka=1
      kb=ka+ncontractions**2
      kc=kb+ncontractions
      kd=kc+ncontractions
      ke=kd+ncontractions**2
      kf=ke+ncontractions**2
      kg=kf+ncontractions**2
      kh=kf+ncontractions
      ki=kf+ncontractions

      if(ki.gt.nreal) stop 'not enough memory for annihilation routine'

      rewind 51

      if(nalpha.eq.nbeta) then
                                call annihilate(nmaxorbs,ncontractions,
     &                                          nalpha,iwkvec,
     &                                          iwkvec(3*numprocs+1),
     &                                          dmaxpc,tol,dmta,focka,
     &                                          coeffa,wkvec(ka),
     &                                          wkvec(kb),wkvec(kc),
     &                                          wkvec(kd),wkvec(ke),
     &                                          wkvec(kf),wkvec(kg),
     &                                          wkvec(kh))

                                do 1002 i=1,ndim
 1002                           dmtb(i)=dmta(i)
                          else
                                call annihilate(nmaxorbs,ncontractions,
     &                                          nalpha,iwkvec,
     &                                          iwkvec(3*numprocs+1),
     &                                          dmaxpc,tol,dmta,focka,
     &                                          coeffa,wkvec(ka),
     &                                          wkvec(kb),wkvec(kc),
     &                                          wkvec(kd),wkvec(ke),
     &                                          wkvec(kf),wkvec(kg),
     &                                          wkvec(kh))

                                call annihilate(nmaxorbs,ncontractions,
     &                                          nbeta,iwkvec,
     &                                          iwkvec(3*numprocs+1),
     &                                          dmaxpc,tol,dmtb,fockb,
     &                                          coeffb,wkvec(ka),
     &                                          wkvec(kb),wkvec(kc),
     &                                          wkvec(kd),wkvec(ke),
     &                                          wkvec(kf),wkvec(kg),
     &                                          wkvec(kh))
                          endif

      ja=1
      jb=ja+ncontractions*ncontractions
      jc=jb+ncontractions*ncontractions
      jd=jc+ncontractions

      call totale(0,ncontractions,ncdfunc,iteration,nptcharges,natoms,
     &            nalpha,nbeta,ncdfunc,energy,diiserror,dmta,dmtb,focke,
     &            coord,charge,coordpc,ptcharge,wkvec,cdfitc,dmaxpc)

      if(nxcfunc.gt.0) then

            do 1003 i=1,ndim
 1003       wkvec(i)=dmta(i)+dmtb(i)

            call xvecmm(nreal-ndim,ninteger,nconts,ncontp,ncontd,nxcs,
     &                  nxcspd,nshels,nshelp,nsheld,mtloca,ilfunc,
     &                  ngaussians,icfunc,ixclfunc,ixccfunc,iwkvec,
     &                  cutoff(3),coord,alpha,coeffs,alphaxc,coefsxc,
     &                  coefpxc,coefdxc,wkvec(1),tveca,wkvec(ndim+1))

                       endif

      nnn=max(nxcfunc,ncontractions)

      ja=1
      jb=ja+3*nnn
      jc=jb+3*nnn
      jd=jc+nnn*nnn

      call gridcall(nreal-jd,nnn,1,iteration+1,nfunctional,natoms,
     &              nmaxatoms,nbsystems,nalpha,nbeta,nconts,ncontp,
     &              ncontd,nxcs,nxcs+nxcspd,nmaxpoints,natomtype,
     &              nptsatom,nshels,nshelp,nsheld,ilfunc,icfunc,
     &              ngaussians,ixclfunc,ixccfunc,isys,iwkvec,
     &              iwkvec(nnn+1),dmixing,cutoff(2),coord,alpha,coeffs,
     &              alphaxc,coefsxc,coefpxc,coefdxc,xcfita,xcfitb,
     &              xcfite,dmta,dmtb,tveca,tvecb,tvece,wkvec(ja),
     &              wkvec(jb),wkvec(jc),focka,fockb,focke,wkvec(jd))

      if(myid.eq.0) then
                          rewind 43
                          call fastrd(43,overlap,ndim)
                    endif

      call mpi_bcast(overlap,ndim,mpi_real8,0,mpi_comm_world,ierr)

      if(myid.eq.0) call cdpops(level,ncontractions,natoms,nmaxatoms,
     &                          nconts,ncontp,ncontd,nbsystems,nshels,
     &                          nshelp,nsheld,icfunc,nfuncatom,isys,
     &                          overlap,dmta,dmtb,wkvec(1),
     &                          wkvec(nbsystems+0*(ncontractions**2)+1),
     &                          wkvec(nbsystems+1*(ncontractions**2)+1),
     &                          wkvec(nbsystems+2*(ncontractions**2)+1),
     &                          wkvec(nbsystems+3*(ncontractions**2)+1))

      do 1004 i=1,ncdfunc
      storecd(i)=zero
 1004 cdstore(i)=cdfitc(i)

      ncds_save=ncds
      ncdspd_save=ncdspd

      do 1005 i=1,ncdfunc
 1005 icdlsave(i)=icdlfunc(i)

      if(myid.eq.0) rewind 18

      do 1006 i=1,nbsystems

      do 1007 j=1,ndim
 1007 wkvec(nbsystems+j)=dmta(j)+dmtb(j)

      do 1008 j=1,natoms
 1008 iwkvec(j)=0

      do 1009 j=1,natoms
      do 1009 k=1,isys(1,i)
 1009 if(dist(coord(1,j),coord(1,isys(k+1,i))).le.cutoff(1)) iwkvec(j)=1

      ncds=0

      do 1010 j=1,ncds_save
      if(iwkvec(icdcfunc(icdlsave(j))).eq.1) then

                             ncds=ncds+1
                             icdlfunc(ncds)=icdlsave(j)

                                             endif
 1010 continue

      ncdspd=0

      do 1011 j=1,ncdspd_save
      if(iwkvec(icdcfunc(icdlsave(ncds_save+j))).eq.1) then

                             ncdspd=ncdspd+1
                             icdlfunc(ncds+ncdspd)=icdlsave(ncds_save+j)

                                                       endif
 1011 continue

      do 1012 j=1,natoms
 1012 iwkvec(j)=0

      do 1013 k=1,level
      do 1013 j=1,isys(1,i)
 1013 if(myid.eq.0) iwkvec((k-1)*(natoms/level)+isys(j+1,i))=1

      call mpi_bcast(iwkvec,natoms,mpi_integer,0,mpi_comm_world,ierr)

      call cdftmm(nreal-nbsystems-ndim,ninteger,nconts,ncontp,ncontd,
     &            ncds,ncdspd,nshels,nshelp,nsheld,mtloca,ilfunc,
     &            ngaussians,icfunc,icdlfunc,icdcfunc,iwkvec,
     &            iwkvec(nmaxatoms+1),coord,alpha,coeffs,alphacd,
     &            coefscd,coefpcd,coefdcd,tvector,wkvec(nbsystems+1),
     &            wkvec(nbsystems+ndim+1))

      nsize=ncds+10*ncdspd

      if(myid.eq.0) call docdfit(nsize,ncds,ncdspd,icdlfunc,wkvec(i),
     &                           alphacd,coefscd,coefdcd,cdfitc,tvector,
     &                           wkvec(nbsystems+1),
     &                           wkvec(nbsystems+nsize*nsize+1))

      do 1014 j=1,ncdfunc
 1014 tvector(j)=cdfitc(j)

      do 1015 j=1,natoms
 1015 iwkvec(j)=0

      do 1016 j=1,natoms
      do 1016 k=1,isys(1,i)
 1016 if(dist(coord(1,j),coord(1,isys(k+1,i))).le.cutoff(1)) iwkvec(j)=1

      ncds=0

      do 1017 j=1,ncds_save
      if(iwkvec(icdcfunc(icdlsave(j))).eq.1) then

                                             ncds=ncds+1
                                             cdfitc(j)=tvector(ncds)

                                             else

                                             cdfitc(j)=zero

                                             endif
 1017 continue

      ncdspd=0

      do 1018 j=1,ncdspd_save
      if(iwkvec(icdcfunc(icdlsave(ncds_save+j))).eq.1) then

            ncdspd=ncdspd+1
            cdfitc(ncds_save+10*(j-1)+01)=tvector(ncds+10*(ncdspd-1)+01)
            cdfitc(ncds_save+10*(j-1)+02)=tvector(ncds+10*(ncdspd-1)+02)
            cdfitc(ncds_save+10*(j-1)+03)=tvector(ncds+10*(ncdspd-1)+03)
            cdfitc(ncds_save+10*(j-1)+04)=tvector(ncds+10*(ncdspd-1)+04)
            cdfitc(ncds_save+10*(j-1)+05)=tvector(ncds+10*(ncdspd-1)+05)
            cdfitc(ncds_save+10*(j-1)+06)=tvector(ncds+10*(ncdspd-1)+06)
            cdfitc(ncds_save+10*(j-1)+07)=tvector(ncds+10*(ncdspd-1)+07)
            cdfitc(ncds_save+10*(j-1)+08)=tvector(ncds+10*(ncdspd-1)+08)
            cdfitc(ncds_save+10*(j-1)+09)=tvector(ncds+10*(ncdspd-1)+09)
            cdfitc(ncds_save+10*(j-1)+10)=tvector(ncds+10*(ncdspd-1)+10)

                                                       else

            cdfitc(ncds_save+10*(j-1)+01)=zero
            cdfitc(ncds_save+10*(j-1)+02)=zero
            cdfitc(ncds_save+10*(j-1)+03)=zero
            cdfitc(ncds_save+10*(j-1)+04)=zero
            cdfitc(ncds_save+10*(j-1)+05)=zero
            cdfitc(ncds_save+10*(j-1)+06)=zero
            cdfitc(ncds_save+10*(j-1)+07)=zero
            cdfitc(ncds_save+10*(j-1)+08)=zero
            cdfitc(ncds_save+10*(j-1)+09)=zero
            cdfitc(ncds_save+10*(j-1)+10)=zero

                                             endif
 1018 continue

      do 1019 j=1,ncdfunc
 1019 storecd(j)=storecd(j)+cdfitc(j)

      ncds=ncds_save
      ncdspd=ncdspd_save

      do 1020 j=1,ncdfunc
 1020 icdlfunc(j)=icdlsave(j)

 1006 continue

      do 1021 i=1,ncdfunc
 1021 cdfitc(i)=dmixing*storecd(i)+(one-dmixing)*cdstore(i)

      return
      end
