      subroutine delphidr(nreal,ninteger,natoms,nmaxorbs,nmaxaux,
     &                    nmaxpts,nmaxesp,mesp,nmaxiterations,idiis,
     &                    ngridtype,nfunctional,iptsflag,ncenters,
     &                    nptcharges,nbsystems,level,nalpha,nbeta,
     &                    nconts,ncontp,ncontd,ncds,ncdspd,nxcs,nxcspd,
     &                    ixcfit,natomtype,nptsatom,nshels,nshelp,
     &                    nsheld,ilfunc,icfunc,ngaussians,mtloca,
     &                    nfuncatom,icdlfunc,icdcfunc,ixclfunc,ixccfunc,
     &                    isys,iwkvec,cutoff,epsilon,dconvergence,
     &                    dmixing,elevelshift,coord,charge,
     &                    coordptcharge,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,xcfitsave,
     &                    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)

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

      dimension cutoff(*),coord(3,*),charge(*),coordptcharge(3,*),
     &          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(*),
     &          xcfitsave(*),wkvec(*)

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

      ncontractions=nconts+3*ncontp+6*ncontd

      ncdtotal=ncds+10*ncdspd
      nxctotal=nxcs+10*nxcspd

      savedmixing=dmixing

      ii=0
      ndelphi=0

      do 1001 i=1,nmaxiterations

      ii=ii+1

      if(ndelphi.eq.5) goto 1001

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

      call scfcycle(ii,nreal-jc,ninteger-nmaxaux,natoms,nmaxorbs,
     &              nmaxaux,nmaxpts,idiis,nfunctional,ngridtype,
     &              ncenters,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,iwkvec(1),
     &              iwkvec(nmaxaux+1),energy,elevelshift,dmixing,
     &              diiserror,dmaxpc,cutoff,coord,coordptcharge,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,
     &              wkvec(ja),wkvec(jb),wkvec(jc))

      factor=max(dconvergence,dconvergence*(10**(3-ndelphi)))

      if((ii.le.2).or.(dmaxpc.gt.factor)) goto 1001

      ii=0

      write(6,1002)
 1002 format(//)

      ja=1
      jb=ja+3*mesp
      jc=jb+4*natoms*natoms
      jd=jc+2*natoms
      je=jd+natoms
      jf=je+3*natoms
      jg=jf+3*natoms
      jh=jg+3*nmaxesp

      ka=1
      kb=ka+mesp
      kc=kb+mesp
      kd=kb+natoms*2
      ke=kb+natoms*natoms

      call mespgen(nreal-jh,ninteger-ke,iptsflag,natoms,mesp,nmaxesp,
     &             nmaxaux,level,ncenters,nbsystems,0,0,nalpha,nbeta,
     &             nconts,ncontp,ncontd,ncds,ncdspd,nshels,nshelp,
     &             nsheld,ilfunc,ngaussians,icfunc,mtloca,icdlfunc,
     &             icdcfunc,natomtype,iwkvec(kc),iwkvec(kd),isys,
     &             iwkvec(ka),iwkvec(kb),iwkvec(ke),coord,charge,dmta,
     &             dmtb,alpha,coeffs,alphacd,coefscd,coefpcd,coefdcd,
     &             cdfitc,wkvec(ja),wkvec(jb),wkvec(jc),wkvec(jd),
     &             wkvec(je),wkvec(jf),wkvec(jg),wkvec(jh))

      ndelphi=ndelphi+1

      ngrid=65
      ngridby2=33

      do 1003 j=1,ncenters
 1003 wkvec(j)=wkvec(jc+j-1)

      ja=1
      jb=ja+ncenters
      jc=jb+ngrid*ngrid*ngrid
      jd=jc+ngrid*ngrid*ngrid
      je=jd+ngrid*ngrid*ngrid
      jf=je+ngrid*ngrid*ngrid*6
      jg=jf+ngrid*ngrid*ngrid
      jh=jg+ncenters*3

      ji=jh+ncenters

      if(ji.gt.nreal) stop 'delphi driver real overflow'

      if((ngrid**3).gt.ninteger) stop 'delphi driver integer overflow'

      if(ndelphi.eq.5) iswitch=1

      call delphi(iswitch,ndelphi,ncenters,nptcharges,ngrid,ngridby2,
     &            natomtype,iwkvec,epsilon,coord,wkvec(ja),
     &            coordptcharge,ptcharge,wkvec(jb),wkvec(jc),wkvec(jd),
     &            wkvec(je),wkvec(jf),wkvec(jg),wkvec(jh))

      ja=1
      jb=ja+(ncenters+nptcharges)*3
      jc=jb+(ncenters+nptcharges)

      call coremm(nreal-jc,ninteger,ncenters,nptcharges,nconts,ncontp,
     &            ncontd,nshels,nshelp,nsheld,mtloca,ilfunc,ngaussians,
     &            icfunc,iwkvec,cutoff(3),wkvec(ja),coord,coordptcharge,
     &            wkvec(jb),charge,ptcharge,alpha,coeffs,overlap,core,
     &            wkvec(jc))

 1001 continue

      dmixing=savedmixing

      if(dmaxpc.le.dconvergence) then
                                       write(6,1004)
                                 else
                                       write(6,1005)
                                 endif

 1004 format(/,' scf converged!!! :-) :-) :-) ')
 1005 format(/,' scf terminated unsuccessfully??? :-( :-( :-( ')

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

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

      call gridcall(nreal-jd,nxctotal,0,-1,nfunctional,ncenters,natoms,
     &              nbsystems,nalpha,nbeta,nconts,ncontp,ncontd,nxcs,
     &              nxcs+nxcspd,nmaxpts,natomtype,nptsatom,nshels,
     &              nshelp,nsheld,ilfunc,icfunc,ngaussians,ixclfunc,
     &              ixccfunc,isys,iwkvec,iwkvec(nxctotal+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 1006 i=1,nxctotal
 1006 xcfitsave(i)=xcfite(i)

      do 1007 i=1,nxctotal
 1007 xcfite(i)=zero

      call numxcc(nreal,ninteger,natoms,nmaxpts,ncenters,nconts,ncontp,
     &            ncontd,ncds,ncdspd,nxcs,nxcspd,nbsystems,level,iwkvec,
     &            mtloca,ilfunc,ngaussians,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 1008 i=1,nxctotal
 1008 xcfite(i)=xcfitsave(i)

      call totale(1,ncontractions,ncdtotal,-1,nptcharges,ncenters,
     &            nalpha,nbeta,nmaxaux,energy,diiserror,dmta,dmtb,focke,
     &            coord,charge,coordptcharge,ptcharge,wkvec,cdfitc,
     &            dmaxpc)

      return
      end
