      subroutine scfrun(nreal,ninteger,nmaxatoms,nmaxorbs,nmaxaux,
     &                  nmaxpoints,nfunctional,maxiterations,idiis,
     &                  ngridtype,natoms,nptcharges,nbsystems,level,
     &                  nalpha,nbeta,nconts,ncontp,ncontd,ncds,ncdspd,
     &                  nxcs,nxcspd,ixcfit,natomtype,nptsatom,nshels,
     &                  nshelp,nsheld,ilfunc,icfunc,mtloca,ngaussian,
     &                  nfuncatom,icdlfunc,icdcfunc,ixclfunc,ixccfunc,
     &                  isys,iwkvec,converged,dmixing,eshift,diiserror,
     &                  cutoff,coord,charge,coordpc,ptcharge,core,focka,
     &                  fockb,focke,dmta,dmtb,coeffa,coeffb,alpha,
     &                  coeffs,alphacd,coefscd,coefpcd,coefdcd,alphaxc,
     &                  coefsxc,coefpxc,coefdxc,cdfitc,xcfita,xcfitb,
     &                  xcfite,xcfite_save,tvector,tveca,tvecb,tvece,
     &                  overlapinv,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(*),
     &          ilfunc(*),icfunc(*),mtloca(*),ngaussian(*),nfuncatom(*),
     &          icdlfunc(*),icdcfunc(*),ixclfunc(*),ixccfunc(*),
     &          isys(nmaxatoms,*),iwkvec(*)

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

      data zero,one,dlarge/0.0d0,1.0d0,1.0d+25/

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

      call mpi_bcast(natoms,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(nptcharges,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(nbsystems,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(level,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(nalpha,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(nbeta,1,mpi_integer,0,mpi_comm_world,ierr)

      call mpi_bcast(nconts,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(ncontp,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(ncontd,1,mpi_integer,0,mpi_comm_world,ierr)

      call mpi_bcast(ncds,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(ncdspd,1,mpi_integer,0,mpi_comm_world,ierr)

      call mpi_bcast(nxcs,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(nxcspd,1,mpi_integer,0,mpi_comm_world,ierr)

      ncontractions=nconts+3*ncontp+6*ncontd

      ndim=ncontractions*(ncontractions+1)/2

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

      call mpi_bcast(idiis,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(ngridtype,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(nfunctional,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(maxiterations,1,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(ixcfit,1,mpi_integer,0,mpi_comm_world,ierr)

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

      call mpi_bcast(nshels,nconts,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(nshelp,ncontp,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(nsheld,ncontd,mpi_integer,0,mpi_comm_world,ierr)

      call mpi_bcast(mtloca,ncontractions,
     &               mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(ilfunc,ncontractions,
     &               mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(icfunc,ncontractions,
     &               mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(ngaussian,ncontractions,
     &               mpi_integer,0,mpi_comm_world,ierr)

      call mpi_bcast(icdlfunc,ncdfunc,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(icdcfunc,ncdfunc,mpi_integer,0,mpi_comm_world,ierr)

      call mpi_bcast(ixclfunc,nxcfunc,mpi_integer,0,mpi_comm_world,ierr)
      call mpi_bcast(ixccfunc,nxcfunc,mpi_integer,0,mpi_comm_world,ierr)

      call mpi_bcast(isys,nmaxatoms*nmaxatoms,
     &               mpi_integer,0,mpi_comm_world,ierr)

      call mpi_bcast(converged,1,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(dmixing,1,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(cutoff(1),3,mpi_real8,0,mpi_comm_world,ierr)

      call mpi_bcast(coord,3*natoms,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(charge,natoms,mpi_real8,0,mpi_comm_world,ierr)

      call mpi_bcast(coordpc,3*nptcharges,
     &               mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(ptcharge,nptcharges,
     &               mpi_real8,0,mpi_comm_world,ierr)

      call mpi_bcast(alpha,10000,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(coeffs,10000,mpi_real8,0,mpi_comm_world,ierr)

      call mpi_bcast(cdfitc,ncdfunc,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(xcoefa,nxcfunc,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(xcoefb,nxcfunc,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(xcoefe,nxcfunc,mpi_real8,0,mpi_comm_world,ierr)

      call mpi_bcast(alphacd,ncdfunc,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(coefscd,ncdfunc,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(coefpcd,ncdfunc,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(coefdcd,ncdfunc,mpi_real8,0,mpi_comm_world,ierr)

      call mpi_bcast(alphaxc,nxcfunc,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(coefsxc,nxcfunc,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(coefpxc,nxcfunc,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(coefdxc,nxcfunc,mpi_real8,0,mpi_comm_world,ierr)

      e1=zero
      e2=zero

      dmixing_save=dmixing
      eshift_save=eshift
      idiis_save=idiis

      deltanrg=dlarge
      dmaxpc=dlarge

      ihit=0

      do 1001 i=1,maxiterations

      if(deltanrg.le.converged) then
                                      ihit=ihit+1
                                else
                                      ihit=0
                                endif

      if((ihit.ge.2).and.(dmaxpc.le.converged)) goto 1001

      ja=1
      jb=ja+nmaxaux*3
      jc=jb+nmaxaux*3

      call scfcycle(i,nreal-jc,ninteger-nmaxaux,nmaxatoms,nmaxorbs,
     &              nmaxaux,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,ngaussian,nfuncatom,
     &              icdlfunc,icdcfunc,ixclfunc,ixccfunc,isys,iwkvec(1),
     &              iwkvec(nmaxaux+1),energy,eshift,dmixing,diiserror,
     &              dmaxpc,cutoff,coord,coordpc,charge,ptcharge,core,
     &              focka,fockb,focke,dmta,dmtb,coeffa,coeffb,
     &              overlapinv,alpha,coeffs,alphacd,coefscd,coefpcd,
     &              coefdcd,alphaxc,coefsxc,coefpxc,coefdxc,cdfitc,
     &              xcfita,xcfitb,xcfite,tvector,tveca,tvecb,tvece,
     &              wkvec(ja),wkvec(jb),wkvec(jc))

      e2=e1
      e1=energy

      deltanrg=abs(e2-e1)

      if(dmaxpc.le.converged) then
                                    dmixing=zero
                                    idiis=0
                              endif

      call mpi_bcast(deltanrg,1,mpi_real8,0,mpi_comm_world,ierr)

 1001 continue

      dmixing=dmixing_save
      idiis=idiis_save

      if(deltanrg.le.converged) then
                                      if(myid.eq.0) write(6,1002) 
                                else
                                      if(myid.eq.0) write(6,1003)
                                endif

 1002 format(/,' scf converged!!! :-) :-) :-) ')
 1003 format(/,' scf terminated unsuccessfully??? :-( :-( :-( ')

      if(nalpha.eq.nbeta) then
                                do 1004 j=1,nbeta
                                do 1004 i=1,ncontractions
 1004                           coeffb(i,j)=coeffa(i,j)
                          endif

      call gengrid(nreal,nmaxpoints,ngridtype,0,natoms,nconts,ncontp,
     &             ncontd,ilfunc,icfunc,ngaussian,nshels,nshelp,nsheld,
     &             natomtype,nptsatom,iwkvec,coord,alpha,coeffs,wkvec)

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

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

      call gridcall(nreal-jd,nxcfunc,0,-1,nfunctional,natoms,nmaxatoms,
     &              nbsystems,nalpha,nbeta,nconts,ncontp,ncontd,nxcs,
     &              nxcs+nxcspd,nmaxpoints,natomtype,nptsatom,nshels,
     &              nshelp,nsheld,ilfunc,icfunc,ngaussian,ixclfunc,
     &              ixccfunc,isys,iwkvec,iwkvec(nxcfunc+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))

      do 1005 i=1,nxcfunc
 1005 xcfite_save(i)=xcfite(i)

      do 1006 i=1,nxcfunc
 1006 xcfite(i)=zero

      call numxcc(nreal,ninteger,nmaxatoms,nmaxpoints,natoms,nconts,
     &            ncontp,ncontd,ncds,ncdspd,nxcs,nxcspd,nbsystems,level,
     &            iwkvec,mtloca,ilfunc,ngaussian,icfunc,nshels,nshelp,
     &            nsheld,natomtype,nptsatom,icdlfunc,icdcfunc,ixclfunc,
     &            ixccfunc,isys,cutoff(3),wkvec,coord,core,focka,fockb,
     &            focke,alpha,coeffs,alphacd,coefscd,coefpcd,coefdcd,
     &            alphaxc,coefsxc,coefpxc,coefdxc,cdfitc,xcfita,xcfitb,
     &            xcfite)

      do 1007 i=1,nxcfunc
 1007 xcfite(i)=xcfite_save(i)

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

      return
      end
