      subroutine startscf(nreal,ninteger,nmaxatoms,nmaxorbs,nmaxaux,
     &                    nmaxpoints,natoms,nptcharges,nconts,ncontp,
     &                    ncontd,ncds,ncdspd,nxcs,nxcspd,level,
     &                    igridtype,nbsystems,isys,icdlsave,iwkvec,
     &                    nshels,nshelp,nsheld,ilfunc,icfunc,ngaussians,
     &                    matloca,icdlfunc,icdcfunc,ixclfunc,ixccfunc,
     &                    natomtype,nptsatom,energy,rhochange,diiserror,
     &                    cutoff,wkvec,core,focka,fockb,focke,overlap,
     &                    coeffa,coeffb,coord,coordptcharges,charge,
     &                    ptcharge,alpha,coeff,alphacd,coefscd,coefpcd,
     &                    coefdcd,alphaxc,coefsxc,coefpxc,coefdxc,
     &                    cdcoef,xcoefa,xcoefb,xcoefe,dmta,dmtb,tveca,
     &                    tvecb,tvece)

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 isys(nmaxatoms,*),icdlsave(*),iwkvec(*),nshels(*),
     &          nshelp(*),nsheld(*),ilfunc(*),icfunc(*),ngaussians(*),
     &          matloca(*),icdlfunc(*),icdcfunc(*),ixclfunc(*),
     &          ixccfunc(*),natomtype(*),nptsatom(*)

      dimension cutoff(*),wkvec(*),core(*),focka(*),fockb(*),focke(*),
     &          overlap(nmaxorbs,*),coeffa(nmaxorbs,*),
     &          coeffb(nmaxorbs,*),coord(3,*),coordptcharges(3,*),
     &          charge(*),ptcharge(*),alpha(*),coeff(*),alphacd(*),
     &          coefscd(*),coefpcd(*),coefdcd(*),alphaxc(*),coefsxc(*),
     &          coefpxc(*),coefdxc(*),cdcoef(*),xcoefa(*),xcoefb(*),
     &          xcoefe(*),dmta(*),dmtb(*),tveca(*),tvecb(*),tvece(*)

      data zero,one,bignumber/0.0d0,1.0d0,1.0d+35/

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

      if(myid.eq.0) then
                          rewind 43
                          rewind 45
                          rewind 66
                    endif

      ncontractions=nconts+3*ncontp+6*ncontd

      ndim=ncontractions*(ncontractions+1)/2

      energy=zero
      rhochange=zero
      diiserror=zero

      call gengrid(nreal,nmaxpoints,igridtype,1,natoms,nconts,ncontp,
     &             ncontd,ilfunc,icfunc,ngaussians,nshels,nshelp,nsheld,
     &             natomtype,nptsatom,iwkvec,coord,alpha,coeff,wkvec)

      ja=1
      jb=ja+(natoms+nptcharges)*3
      jc=jb+(natoms+nptcharges)
      jd=jc+ndim

      call coremm(nreal-jd-ndim,ninteger,natoms,nptcharges,nconts,
     &            ncontp,ncontd,nshels,nshelp,nsheld,matloca,ilfunc,
     &            ngaussians,icfunc,iwkvec,cutoff(3),wkvec(ja),coord,
     &            coordptcharges,wkvec(jb),charge,ptcharge,alpha,coeff,
     &            wkvec(jc),core,wkvec(jd+1))

      if(myid.eq.0) call fastwr(43,wkvec(jc),ndim)

      do 1001 j=1,ncontractions
      do 1001 k=1,ncontractions
      jk=((max(j,k)*(max(j,k)-1))/2)+min(j,k)-1
 1001 overlap(j,k)=wkvec(jc+jk)

      if(myid.eq.0) then
                          do 1002 i=1,ncontractions
 1002                     call fastwr(45,overlap(1,i),ncontractions)
                    endif

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

      call s_annihilate(nmaxorbs,ncontractions,iwkvec,
     &                  iwkvec(numprocs+1),overlap,coeffa,coeffb,
     &                  wkvec(ka),wkvec(kb),wkvec(kc),wkvec(kd),
     &                  wkvec(ke),wkvec(kf),wkvec(kg))

      nsize=ncds+10*ncdspd

      call cdcalls(nreal-nsize*nsize,ncds,ncdspd,nsize,nmaxatoms,natoms,
     &             nbsystems,isys,icdlsave,icdlfunc,icdcfunc,iwkvec,
     &             cutoff(1),alphacd,coefscd,coefpcd,coefdcd,wkvec,
     &             coord,wkvec(nsize*nsize+1))

      return
      end
