      subroutine gridcall(nreal,nxcdim,iswitch,iteration,nfunctional,
     &                    ncenters,natoms,nbsystems,nalpha,nbeta,nconts,
     &                    ncontp,ncontd,nxcs,nxcspd,nmaxpts,natomtype,
     &                    nptsatom,nshels,nshelp,nsheld,ilfunc,icfunc,
     &                    ngauss,ixclfunc,ixccfunc,isys,ixclsave,iwkvec,
     &                    xim,cutoff,coord,alpha,coeffs,alphaxc,coefsxc,
     &                    coefpxc,coefdxc,xcfita,xcfitb,xcfite,amat,
     &                    bmat,tveca,tvecb,tvece,storexc,xcstore,
     &                    overlap,focka,fockb,focke,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(*),ngauss(*),ixclfunc(*),ixccfunc(*),
     &          isys(natoms,*),ixclsave(*),iwkvec(*)

      dimension coord(3,*),amat(*),bmat(*),alpha(*),coeffs(*),
     &          alphaxc(*),coefsxc(*),coefpxc(*),coefdxc(*),xcfita(*),
     &          xcfitb(*),xcfite(*),tveca(*),tvecb(*),tvece(*),
     &          storexc(*),xcstore(*),overlap(nxcdim,*),focka(*),
     &          fockb(*),focke(*),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)

      ncntrt=nconts+3*ncontp+6*ncontd
      ndim=ncntrt*(ncntrt+1)/2

      nxcfit=nxcs+10*(nxcspd-nxcs)

      if(myid.eq.0) rewind 19

casa
c     if((iteration.gt.1).and.(myid.eq.0)) then

c                              do 1001 i=1,nxcfit
c1001                          call fastrd(19,overlap(1,i),nxcfit)

c                              do 1002 i=1,nxcfit
c                              tvecb(i)=zero
c                              do 1002 j=1,nxcfit
c1002                          tvecb(i)=tvecb(i)+overlap(j,i)*tveca(j)

c                              rewind 19

c                                          endif

c     call mpi_bcast(tvecb,nxcfit,mpi_real8,0,mpi_comm_world,ierr)
casa

      nbsystems_save=nbsystems

      if(nxcfit.eq.0) nbsystems=1

      do 1003 i=1,ndim
      focka(i)=zero
      fockb(i)=zero
 1003 if(nxcfit.eq.0) focke(i)=zero

      call mpi_bcast(amat,ndim,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(bmat,ndim,mpi_real8,0,mpi_comm_world,ierr)

casa
c     do 1004 i=1,nxcfit
c     xcfita(i)=zero
c     xcfitb(i)=zero
c     storexc(2*nxcfit+i)=zero
c1004 xcstore(2*nxcfit+i)=xcfite(i)
casa

      nxcs_save=nxcs
      nxcspd_save=nxcspd-nxcs

      do 1005 i=1,nxcfit
 1005 ixclsave(i)=ixclfunc(i)

      do 1006 i=1,nbsystems

      if(iswitch.eq.0) goto 1007

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

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

      nxcs=0

      do 1010 j=1,nxcs_save
      if(iwkvec(ixccfunc(ixclsave(j))).eq.1) then

                             nxcs=nxcs+1
                             ixclfunc(nxcs)=ixclsave(j)

                                             endif
 1010 continue

      nxcspd=0

      do 1011 j=1,nxcspd_save
      if(iwkvec(ixccfunc(ixclsave(nxcs_save+j))).eq.1) then

                             nxcspd=nxcspd+1
                             ixclfunc(nxcs+nxcspd)=ixclsave(nxcs_save+j)

                                                       endif
 1011 continue

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

      do 1013 j=1,isys(1,i)
 1013 iwkvec(isys(j+1,i))=1

 1007 continue

      nxcfunc=nxcs+10*nxcspd

      nxcspd=nxcspd+nxcs

      nstart=min(512,32*numprocs)

      npoints=nstart+16

 1014 npoints=npoints-16

      if(npoints.le.0) stop 'real overflow for scf grid work'

      nfunctions=max(ncntrt,nxcfunc)
      norbitals=max(nalpha,nbeta)
      nwork=max((nfunctions*30*2+nfunctions*3),(npoints*121),
     &          (((nxcfunc*(nxcfunc+1))/2)+npoints*7))

      ka=1
      kb=ka+npoints*nfunctions
      kc=kb+npoints*norbitals
      kd=kc+npoints*norbitals
      ke=kd+npoints*3
      kf=ke+nmaxpts*2
      kg=kf+nmaxpts
      kh=kg+npoints
      ki=kh+npoints*3
      kj=ki+npoints*3
      kk=kj+npoints*3
      kl=kk+npoints*3
      km=kl+npoints*3
      kn=km+npoints*3
      ko=kn+npoints*3
      kp=ko+npoints*3
      kq=kp+npoints*3
      kr=kq+nmaxpts
      ks=kr+nmaxpts
      kt=ks+nmaxpts
      ku=kt+nmaxpts
      kv=ku+nwork
      kw=kv+npoints*nfunctions
      kx=kw+npoints*nfunctions
      ky=kx+npoints*ncntrt
      kz=ky+npoints*ncntrt
      la=kz+npoints*ncntrt
      lb=la+npoints*ncntrt
      lc=lb+npoints*ncntrt
      ld=lc+npoints*ncntrt
      le=ld+npoints*ncntrt
      lf=le+npoints*ncntrt
      lg=lf+npoints*ncntrt
      lh=lg+npoints*ncntrt
      li=lh+npoints*ncntrt
      lj=li+ncntrt
      lk=lj+ncntrt+nfunctions
      ll=lk+ncntrt+nfunctions
      lm=ll+ncenters
      ln=lm+natoms*natoms

      if((ln.gt.nreal).and.(npoints.gt.16)) goto 1014

      call gridwork(i,iswitch,iteration,nfunctional,nxcdim,npoints,
     &              nfunctions,norbitals,ncntrt,ncenters,natoms,
     &              nbsystems,nalpha,nbeta,nconts,ncontp,ncontd,nxcfunc,
     &              nxcs,nxcspd,natomtype,nptsatom,nshels,nshelp,nsheld,
     &              ilfunc,icfunc,ngauss,ixclfunc,ixccfunc,iwkvec,
     &              iwkvec(ncenters+1),coord,alpha,coeffs,alphaxc,
     &              coefsxc,coefpxc,coefdxc,tveca,tvecb,tvece,overlap,
     &              amat,bmat,wkvec(ka),wkvec(kb),wkvec(kc),wkvec(kd),
     &              wkvec(ke),wkvec(kf),wkvec(kg),wkvec(kh),wkvec(ki),
     &              wkvec(kj),wkvec(kk),wkvec(kl),wkvec(km),wkvec(kn),
     &              wkvec(ko),wkvec(kp),wkvec(kq),wkvec(kr),wkvec(ks),
     &              wkvec(kt),wkvec(ku),wkvec(kv),wkvec(kw),wkvec(kx),
     &              wkvec(ky),wkvec(kz),wkvec(la),wkvec(lb),wkvec(lc),
     &              wkvec(ld),wkvec(le),wkvec(lf),wkvec(lg),wkvec(lh),
     &              wkvec(li),wkvec(lj),wkvec(lk),wkvec(ll),wkvec(lm),
     &              focka,fockb,focke)

      if(iswitch.eq.0) goto 1015

casa
c     ja=1
c     jb=ja+nxcfunc*nxcfunc
c     jc=jb+nxcfunc*nxcfunc
c     jd=jc+nxcfunc*nxcfunc
c     je=jd+nxcfunc*nxcfunc
c     jf=je+nxcfunc
c     jg=jf+nxcfunc
c     jh=jg+nxcfunc
c     ji=jh+nxcfunc

c     if(ji.gt.nreal) stop 'matrix inversion required for grid work'

c     if((nxcfunc.gt.1).and.(iteration.eq.1)) then

c             do 1016 j=1,nxcfunc
c             do 1016 k=1,j
c1016         overlap(j,k)=overlap(k,j)

c             jk=0

c             do 1017 k=1,nxcfunc
c             do 1017 j=1,nxcfunc
c             jk=jk+1
c1017         wkvec(jk)=overlap(j,k)

c             call mpi_bcast(wkvec,nxcfunc**2,
c    &                       mpi_real8,0,mpi_comm_world,ierr)

c             call x_annihilate(19,nxcfunc,iwkvec,iwkvec(numprocs+1),
c    &                          wkvec(ja),overlap,wkvec(jb),wkvec(jc),
c    &                          wkvec(jd),wkvec(je),wkvec(jf),wkvec(jg),
c    &                          wkvec(jh))

c             goto 1006

c                                             endif
casa

      if((myid.gt.0).or.(nxcfunc.eq.0)) goto 1006

casa
c     call gtxccf(xcfite,tvece,overlap,nxcfunc,nxcdim)

c     do 1018 j=1,nxcfunc
c1018 tvece(j)=xcfite(j)

c     do 1019 j=1,ncenters
c1019 iwkvec(j)=0

c     do 1020 j=1,ncenters
c     do 1020 k=1,isys(1,i)
c1020 if(dist(coord(1,j),coord(1,isys(k+1,i))).le.cutoff) iwkvec(j)=1

c     nxcs=0

c     do 1021 j=1,nxcs_save
c     if(iwkvec(ixccfunc(ixclsave(j))).eq.1) then
c                                                  nxcs=nxcs+1

c                                                  xcfite(j)=tvece(nxcs)
c                                            else
c                                                  xcfite(j)=zero
c                                            endif
c1021 continue

c     nxcspd=0

c     do 1022 j=1,nxcspd_save
c     if(iwkvec(ixccfunc(ixclsave(nxcs_save+j))).eq.1) then

c             nxcspd=nxcspd+1

c             xcfite(nxcs_save+10*(j-1)+01)=tvece(nxcs+10*(nxcspd-1)+01)
c             xcfite(nxcs_save+10*(j-1)+02)=tvece(nxcs+10*(nxcspd-1)+02)
c             xcfite(nxcs_save+10*(j-1)+03)=tvece(nxcs+10*(nxcspd-1)+03)
c             xcfite(nxcs_save+10*(j-1)+04)=tvece(nxcs+10*(nxcspd-1)+04)
c             xcfite(nxcs_save+10*(j-1)+05)=tvece(nxcs+10*(nxcspd-1)+05)
c             xcfite(nxcs_save+10*(j-1)+06)=tvece(nxcs+10*(nxcspd-1)+06)
c             xcfite(nxcs_save+10*(j-1)+07)=tvece(nxcs+10*(nxcspd-1)+07)
c             xcfite(nxcs_save+10*(j-1)+08)=tvece(nxcs+10*(nxcspd-1)+08)
c             xcfite(nxcs_save+10*(j-1)+09)=tvece(nxcs+10*(nxcspd-1)+09)
c             xcfite(nxcs_save+10*(j-1)+10)=tvece(nxcs+10*(nxcspd-1)+10)

c                                                      else

c             xcfite(nxcs_save+10*(j-1)+01)=zero
c             xcfite(nxcs_save+10*(j-1)+02)=zero
c             xcfite(nxcs_save+10*(j-1)+03)=zero
c             xcfite(nxcs_save+10*(j-1)+04)=zero
c             xcfite(nxcs_save+10*(j-1)+05)=zero
c             xcfite(nxcs_save+10*(j-1)+06)=zero
c             xcfite(nxcs_save+10*(j-1)+07)=zero
c             xcfite(nxcs_save+10*(j-1)+08)=zero
c             xcfite(nxcs_save+10*(j-1)+09)=zero
c             xcfite(nxcs_save+10*(j-1)+10)=zero

c                                            endif
c1022 continue

c     do 1023 j=1,nxcfit
c1023 storexc(2*nxcfit+j)=storexc(2*nxcfit+j)+xcfite(j)

 1006 continue

casa
c     do 1024 j=1,nxcfit
c1024 ixclfunc(j)=ixclsave(j)

c     nxcs=nxcs_save
c     nxcspd=nxcs_save+nxcspd_save

c     if(iteration.gt.1) then

c        do 1025 i=1,nxcfit
c1025    xcfite(i)=xim*storexc(2*nxcfit+i)+(one-xim)*xcstore(2*nxcfit+i)

c                        endif
casa

 1015 continue

      nbsystems=nbsystems_save

      return
      end
