      subroutine gridwork(isubsystem,iswitch,iteration,nfunctional,
     &                    nauxfunctions,npoints,naos,nmos,ncontractions,
     &                    ncenters,nmaxatoms,nbsystems,nalpha,nbeta,
     &                    nconts,ncontp,ncontd,nxcfunctions,nxcs,nxcspd,
     &                    natomtype,nptsatom,nshels,nshelp,nsheld,
     &                    ilfunc,icfunc,ngaussians,ixclfunc,ixccfunc,
     &                    isys,iwkvec,coord,alpha,coeff,alphaxc,coefsxc,
     &                    coefpxc,coefdxc,extent,tvecb,tvece,auxoverlap,
     &                    amat,bmat,fv,fa,fb,dnsty,xcpot,xcnrg,weight,
     &                    gradx,grady,gradz,hssxx,hssxy,hssxz,hssyy,
     &                    hssyz,hsszz,xpoint,ypoint,zpoint,wpoint,wkvec,
     &                    fvx,fvy,fvz,fvxx,fvxy,fvxz,fvyy,fvyz,fvzz,wf,
     &                    wfx,wfy,wfz,extent1,extent2,extent3,extent4,
     &                    dcc,focka,fockb,focke)

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 istat(mpi_status_size)

      dimension natomtype(*),nptsatom(*),nshels(*),nshelp(*),nsheld(*),
     &          ilfunc(*),icfunc(*),ngaussians(*),ixclfunc(*),
     &          ixccfunc(*),isys(*),iwkvec(*)

      dimension coord(3,*),alpha(*),coeff(*),alphaxc(*),coefsxc(*),
     &          coefpxc(*),coefdxc(*),extent(*),tvecb(*),tvece(*),
     &          auxoverlap(nauxfunctions,*),amat(*),bmat(*),
     &          fv(npoints,*),fa(nmos,*),fb(nmos,*),dnsty(npoints,*),
     &          xcpot(npoints,*),xcnrg(*),weight(*),gradx(npoints,*),
     &          grady(npoints,*),gradz(npoints,*),hssxx(npoints,*),
     &          hssxy(npoints,*),hssxz(npoints,*),hssyy(npoints,*),
     &          hssyz(npoints,*),hsszz(npoints,*),xpoint(*),ypoint(*),
     &          zpoint(*),wpoint(*),wkvec(*),fvx(npoints,*),
     &          fvy(npoints,*),fvz(npoints,*),fvxx(npoints,*),
     &          fvxy(npoints,*),fvxz(npoints,*),fvyy(npoints,*),
     &          fvyz(npoints,*),fvzz(npoints,*),wf(npoints,*),
     &          wfx(npoints,*),wfy(npoints,*),wfz(npoints,*),extent1(*),
     &          extent2(*),extent3(*),extent4(*),dcc(ncenters,*),
     &          focka(*),fockb(*),focke(*)

      data zero,half,one,two,three/0.0d0,0.5d0,1.0d0,2.0d0,3.0d0/
      data small,cutoff/1.0d-14,1.0d-16/
      data tolerance,cutrho2/1.0d-18,1.0d-20/

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

      ixcfit=0
      if(nxcfunctions.gt.0) ixcfit=1

      igrad=0
      if(iswitch.eq.0) igrad=1

      ipert=0
      if(nfunctional.eq.2) ipert=1

      ipart=0
      if((ixcfit.eq.0).and.(igrad.eq.0)) ipart=1

      do 1001 i=1,ncenters
      do 1001 j=1,ncenters
 1001 dcc(j,i)=dist(coord(1,j),coord(1,i))

      call diffuse2(nconts,ncontp,ncontd,ilfunc,icfunc,ngaussians,
     &              nshels,nshelp,nsheld,cutoff,alpha,coeff,coord,
     &              extent1)

      call diffuse3(nxcs,nxcspd,ixclfunc,ixccfunc,cutoff,alphaxc,
     &              coefsxc,coefpxc,coefdxc,coord,extent2)

      call diffuse4(ncenters,natomtype,cutrho2,extent4)

      nsavefunctional=nfunctional

      if((ipert.eq.1).and.(igrad.eq.1)) nfunctional=1
      if((ipert.eq.1).and.(igrad.eq.0)) nfunctional=0

      if(myid.eq.0) then
                          rewind 11
                          rewind 12
                    endif

casa
c     if(iteration.eq.1) then
c                              do 1002 i=1,nxcfunctions
c                              do 1002 j=1,nxcfunctions
c1002                          auxoverlap(j,i)=zero
c                        endif

c     do 1003 i=1,nxcfunctions
c1003 tvece(i)=zero
casa

      do 1004 iatom=1,ncenters

      if(natomtype(iatom).eq.0) goto 1004

      if(myid.eq.0) then
                          call fastrd(11,xpoint,nptsatom(iatom))
                          call fastrd(11,ypoint,nptsatom(iatom))
                          call fastrd(11,zpoint,nptsatom(iatom))
                          call fastrd(11,wpoint,nptsatom(iatom))
                    endif

      do 1005 ia=1,nptsatom(iatom),npoints

      ib=min(ia+npoints-1,nptsatom(iatom))

      ilast=0
      if((iatom.eq.ncenters).and.(ib.eq.nptsatom(iatom))) ilast=1

      radius=sqrt((coord(1,iatom)-xpoint(ib))**2
     &           +(coord(2,iatom)-ypoint(ib))**2
     &           +(coord(3,iatom)-zpoint(ib))**2)

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

      npts=ib-ia+1

      nmaxloop=(npts/numprocs)+1
      if((nmaxloop*numprocs).lt.npts) nmaxloop=nmaxloop+1

      if(myid.eq.0) then
                          ja=ia

                          islave=0

                          do 1006 i=1,npts,nmaxloop

                          nloop=min(nmaxloop,npts-i+1)

                          call mpi_bsend(nloop,1,mpi_integer,
     &                                   islave,77,mpi_comm_world,ierr)
                          call mpi_bsend(ja,1,mpi_integer,
     &                                   islave,77,mpi_comm_world,ierr)

                          call mpi_bsend(xpoint(ja),nloop,mpi_real8,
     &                                   islave,77,mpi_comm_world,ierr)
                          call mpi_bsend(ypoint(ja),nloop,mpi_real8,
     &                                   islave,77,mpi_comm_world,ierr)
                          call mpi_bsend(zpoint(ja),nloop,mpi_real8,
     &                                   islave,77,mpi_comm_world,ierr)
                          call mpi_bsend(wpoint(ja),nloop,mpi_real8,
     &                                   islave,77,mpi_comm_world,ierr)

                          if(isubsystem.gt.1) then

                                   call fastrd(12,xcpot(1,1),nloop)
                                   call fastrd(12,xcpot(1,2),nloop)
                                   call fastrd(12,xcnrg(1),nloop)

                                   call mpi_bsend(xcpot(1,1),nloop,
     &                                            mpi_real8,islave,77,
     &                                            mpi_comm_world,ierr)
                                   call mpi_bsend(xcpot(1,2),nloop,
     &                                            mpi_real8,islave,77,
     &                                            mpi_comm_world,ierr)
                                   call mpi_bsend(xcnrg(1),nloop,
     &                                            mpi_real8,islave,77,
     &                                            mpi_comm_world,ierr)

                                              endif

                          ja=ja+nloop

                          last_slave=islave
                           
 1006                     islave=islave+1

                          do 1007 i=islave,numprocs-1
                          nloop=0
 1007                     call mpi_bsend(nloop,1,mpi_integer,
     &                                   i,77,mpi_comm_world,ierr)
                    endif

      call mpi_recv(npts,1,mpi_integer,0,77,mpi_comm_world,istat,ierr)

      if(npts.eq.0) goto 1005

      call mpi_recv(ja,1,mpi_integer,0,77,mpi_comm_world,istat,ierr)

      call mpi_recv(xpoint(ja),npts,mpi_real8,
     &              0,77,mpi_comm_world,istat,ierr)
      call mpi_recv(ypoint(ja),npts,mpi_real8,
     &              0,77,mpi_comm_world,istat,ierr)
      call mpi_recv(zpoint(ja),npts,mpi_real8,
     &              0,77,mpi_comm_world,istat,ierr)
      call mpi_recv(wpoint(ja),npts,mpi_real8,
     &              0,77,mpi_comm_world,istat,ierr)

casa
c     if(isubsystem.gt.1) then
c                               call mpi_recv(xcpot(1,1),npts,
c    &                                        mpi_real8,0,77,
c    &                                        mpi_comm_world,istat,ierr)
c                               call mpi_recv(xcpot(1,2),npts,
c    &                                        mpi_real8,0,77,
c    &                                        mpi_comm_world,istat,ierr)
c                               call mpi_recv(xcnrg(1),npts,
c    &                                        mpi_real8,0,77,
c    &                                        mpi_comm_world,istat,ierr)
c                               goto 1038
c                         endif
casa

      if(nfunctional.eq.0) then

                    call compg0(iatom,ncontractions,npts,npoints,
     &                          ncenters,nconts,ncontp,ncontd,ilfunc,
     &                          icfunc,ngaussians,nshels,nshelp,nsheld,
     &                          radius,coord,alpha,coeff,xpoint(ja),
     &                          ypoint(ja),zpoint(ja),fv,extent1,
     &                          extent3,dcc,wkvec(00*npoints+1),
     &                          wkvec(01*npoints+1),wkvec(02*npoints+1),
     &                          wkvec(03*npoints+1),wkvec(04*npoints+1),
     &                          wkvec(05*npoints+1),wkvec(06*npoints+1),
     &                          wkvec(07*npoints+1),wkvec(08*npoints+1),
     &                          wkvec(09*npoints+1),wkvec(10*npoints+1),
     &                          wkvec(11*npoints+1))

                           endif

      if((nfunctional.eq.1).and.(ipart.eq.1)) then

                    call compg1(iatom,ncontractions,npts,npoints,
     &                          ncenters,nconts,ncontp,ncontd,ilfunc,
     &                          icfunc,ngaussians,nshels,nshelp,nsheld,
     &                          radius,coord,alpha,coeff,xpoint(ja),
     &                          ypoint(ja),zpoint(ja),fv,fvx,fvy,fvz,
     &                          extent1,extent3,dcc,wkvec(00*npoints+1),
     &                          wkvec(01*npoints+1),wkvec(02*npoints+1),
     &                          wkvec(03*npoints+1),wkvec(04*npoints+1),
     &                          wkvec(05*npoints+1),wkvec(06*npoints+1),
     &                          wkvec(07*npoints+1),wkvec(08*npoints+1),
     &                          wkvec(09*npoints+1),wkvec(10*npoints+1),
     &                          wkvec(11*npoints+1),wkvec(12*npoints+1),
     &                          wkvec(13*npoints+1),wkvec(14*npoints+1),
     &                          wkvec(15*npoints+1),wkvec(16*npoints+1),
     &                          wkvec(17*npoints+1),wkvec(18*npoints+1),
     &                          wkvec(19*npoints+1),wkvec(20*npoints+1),
     &                          wkvec(21*npoints+1),wkvec(22*npoints+1),
     &                          wkvec(23*npoints+1),wkvec(24*npoints+1))

                                              endif

      if((nfunctional.eq.1).and.(ipart.eq.0)) then

                    call compg2(iatom,ncontractions,npts,npoints,
     &                          ncenters,nconts,ncontp,ncontd,ilfunc,
     &                          icfunc,ngaussians,nshels,nshelp,nsheld,
     &                          radius,coord,alpha,coeff,xpoint(ja),
     &                          ypoint(ja),zpoint(ja),fv,fvx,fvy,fvz,
     &                          fvxx,fvxy,fvxz,fvyy,fvyz,fvzz,extent1,
     &                          extent3,dcc,wkvec(00*npoints+1),
     &                          wkvec(01*npoints+1),wkvec(02*npoints+1),
     &                          wkvec(03*npoints+1),wkvec(04*npoints+1),
     &                          wkvec(05*npoints+1),wkvec(06*npoints+1),
     &                          wkvec(07*npoints+1),wkvec(08*npoints+1),
     &                          wkvec(09*npoints+1),wkvec(10*npoints+1),
     &                          wkvec(11*npoints+1),wkvec(12*npoints+1),
     &                          wkvec(13*npoints+1),wkvec(14*npoints+1),
     &                          wkvec(15*npoints+1),wkvec(16*npoints+1),
     &                          wkvec(17*npoints+1),wkvec(18*npoints+1),
     &                          wkvec(19*npoints+1),wkvec(20*npoints+1),
     &                          wkvec(21*npoints+1),wkvec(22*npoints+1),
     &                          wkvec(23*npoints+1),wkvec(24*npoints+1),
     &                          wkvec(25*npoints+1),wkvec(26*npoints+1),
     &                          wkvec(27*npoints+1),wkvec(28*npoints+1),
     &                          wkvec(29*npoints+1),wkvec(30*npoints+1),
     &                          wkvec(31*npoints+1),wkvec(32*npoints+1),
     &                          wkvec(33*npoints+1),wkvec(34*npoints+1),
     &                          wkvec(35*npoints+1),wkvec(36*npoints+1),
     &                          wkvec(37*npoints+1),wkvec(38*npoints+1),
     &                          wkvec(39*npoints+1))

                                              endif

      do 1008 j=1,2
      do 1008 k=1,npts
 1008 dnsty(k,j)=zero

      if(nfunctional.gt.0) then
                                 do 1009 j=1,2
                                 do 1009 k=1,npts
                                 gradx(k,j)=zero
                                 grady(k,j)=zero
 1009                            gradz(k,j)=zero

                                 if(ipart.eq.0) then
                                                      do 1010 j=1,2
                                                      do 1010 k=1,npts
                                                      hssxx(k,j)=zero
                                                      hssxy(k,j)=zero
                                                      hssxz(k,j)=zero
                                                      hssyy(k,j)=zero
                                                      hssyz(k,j)=zero
 1010                                                 hsszz(k,j)=zero
                                                endif
                           endif

      ihits=0

      do 1011 i=1,ncontractions
      if(extent3(i).gt.tolerance) then
                                        ihits=ihits+1
                                        iwkvec(ihits)=i
                                  endif
 1011 continue

      do 1012 i=1,ihits
      ii=iwkvec(i)

      do 1013 k=1,npts
 1013 wf(k,ii)=zero

      if((nfunctional.gt.0).and.(ipart.eq.0)) then
                                                    do 1014 k=1,npts
                                                    wfx(k,ii)=zero
                                                    wfy(k,ii)=zero
 1014                                               wfz(k,ii)=zero
                                              endif

      do 1012 j=1,ihits
      jj=iwkvec(j)

      ijmax=max(ii,jj)
      ijmin=min(ii,jj)
      ij=(ijmax*(ijmax-1)/2)+ijmin

      if(ii.eq.jj) then
                         factor=amat(ij)*two
                   else
                         factor=amat(ij)
                   endif

      do 1015 k=1,npts
 1015 wf(k,ii)=wf(k,ii)+factor*fv(k,jj)

      if((nfunctional.gt.0).and.(ipart.eq.0)) then

                            do 1016 k=1,npts
                            wfx(k,ii)=wfx(k,ii)+factor*fvx(k,jj)
                            wfy(k,ii)=wfy(k,ii)+factor*fvy(k,jj)
 1016                       wfz(k,ii)=wfz(k,ii)+factor*fvz(k,jj)

                                              endif

 1012 continue
 
      do 1017 i=1,ihits
      ii=iwkvec(i)
      do 1017 k=1,npts
 1017 dnsty(k,1)=dnsty(k,1)+fv(k,ii)*wf(k,ii)

      do 1018 k=1,npts
 1018 dnsty(k,1)=dnsty(k,1)*half

      if(nfunctional.gt.0) then

           do 1019 i=1,ihits
           ii=iwkvec(i)
           do 1019 k=1,npts
           gradx(k,1)=gradx(k,1)+fvx(k,ii)*wf(k,ii)
           grady(k,1)=grady(k,1)+fvy(k,ii)*wf(k,ii)
 1019      gradz(k,1)=gradz(k,1)+fvz(k,ii)*wf(k,ii)

                           endif

      if((nfunctional.gt.0).and.(ipart.eq.0)) then

           do 1020 i=1,ihits
           ii=iwkvec(i)
           do 1020 k=1,npts
           hssxx(k,1)=hssxx(k,1)+fvxx(k,ii)*wf(k,ii)+fvx(k,ii)*wfx(k,ii)
           hssxy(k,1)=hssxy(k,1)+fvxy(k,ii)*wf(k,ii)+fvx(k,ii)*wfy(k,ii)
           hssxz(k,1)=hssxz(k,1)+fvxz(k,ii)*wf(k,ii)+fvx(k,ii)*wfz(k,ii)
           hssyy(k,1)=hssyy(k,1)+fvyy(k,ii)*wf(k,ii)+fvy(k,ii)*wfy(k,ii)
           hssyz(k,1)=hssyz(k,1)+fvyz(k,ii)*wf(k,ii)+fvy(k,ii)*wfz(k,ii)
 1020      hsszz(k,1)=hsszz(k,1)+fvzz(k,ii)*wf(k,ii)+fvz(k,ii)*wfz(k,ii)

                                              endif

      if(nalpha.eq.nbeta) then

         do 1021 k=1,npts
 1021    dnsty(k,2)=dnsty(k,1)

         if(nfunctional.gt.0) then
                                    do 1022 k=1,npts
                                    gradx(k,2)=gradx(k,1)
                                    grady(k,2)=grady(k,1)
 1022                               gradz(k,2)=gradz(k,1)

                                    if(ipart.eq.0) then

                                                  do 1023 k=1,npts
                                                  hssxx(k,2)=hssxx(k,1)
                                                  hssxy(k,2)=hssxy(k,1)
                                                  hssxz(k,2)=hssxz(k,1)
                                                  hssyy(k,2)=hssyy(k,1)
                                                  hssyz(k,2)=hssyz(k,1)
 1023                                             hsszz(k,2)=hsszz(k,1)

                                                   endif
                              endif

         goto 1033

                          endif

      do 1024 i=1,ihits
      ii=iwkvec(i)

      do 1025 k=1,npts
 1025 wf(k,ii)=zero

      if((nfunctional.gt.0).and.(ipart.eq.0)) then
                                                    do 1026 k=1,npts
                                                    wfx(k,ii)=zero
                                                    wfy(k,ii)=zero
 1026                                               wfz(k,ii)=zero
                                              endif

      do 1024 j=1,ihits
      jj=iwkvec(j)

      ijmax=max(ii,jj)
      ijmin=min(ii,jj)
      ij=(ijmax*(ijmax-1)/2)+ijmin

      if(ii.eq.jj) then
                         factor=bmat(ij)*two
                   else
                         factor=bmat(ij)
                   endif

      do 1027 k=1,npts
 1027 wf(k,ii)=wf(k,ii)+factor*fv(k,jj)

      if((nfunctional.gt.0).and.(ipart.eq.0)) then

                            do 1028 k=1,npts
                            wfx(k,ii)=wfx(k,ii)+factor*fvx(k,jj)
                            wfy(k,ii)=wfy(k,ii)+factor*fvy(k,jj)
 1028                       wfz(k,ii)=wfz(k,ii)+factor*fvz(k,jj)

                                              endif

 1024 continue

      do 1029 i=1,ihits
      ii=iwkvec(i)
      do 1029 k=1,npts
 1029 dnsty(k,2)=dnsty(k,2)+fv(k,ii)*wf(k,ii)

      do 1030 k=1,npts
 1030 dnsty(k,2)=dnsty(k,2)*half

      if(nfunctional.gt.0) then

           do 1031 i=1,ihits
           ii=iwkvec(i)
           do 1031 k=1,npts
           gradx(k,2)=gradx(k,2)+fvx(k,ii)*wf(k,ii)
           grady(k,2)=grady(k,2)+fvy(k,ii)*wf(k,ii)
 1031      gradz(k,2)=gradz(k,2)+fvz(k,ii)*wf(k,ii)

                           endif

      if((nfunctional.gt.0).and.(ipart.eq.0)) then

           do 1032 i=1,ihits
           ii=iwkvec(i)
           do 1032 k=1,npts
           hssxx(k,2)=hssxx(k,2)+fvxx(k,ii)*wf(k,ii)+fvx(k,ii)*wfx(k,ii)
           hssxy(k,2)=hssxy(k,2)+fvxy(k,ii)*wf(k,ii)+fvx(k,ii)*wfy(k,ii)
           hssxz(k,2)=hssxz(k,2)+fvxz(k,ii)*wf(k,ii)+fvx(k,ii)*wfz(k,ii)
           hssyy(k,2)=hssyy(k,2)+fvyy(k,ii)*wf(k,ii)+fvy(k,ii)*wfy(k,ii)
           hssyz(k,2)=hssyz(k,2)+fvyz(k,ii)*wf(k,ii)+fvy(k,ii)*wfz(k,ii)
 1032      hsszz(k,2)=hsszz(k,2)+fvzz(k,ii)*wf(k,ii)+fvz(k,ii)*wfz(k,ii)

                                              endif

 1033 continue

      call addvecs(npts,dnsty(1,3),dnsty(1,1),dnsty(1,2))

      if(nfunctional.gt.0) then

                     call addvecs(npts,gradx(1,3),gradx(1,1),gradx(1,2))
                     call addvecs(npts,grady(1,3),grady(1,1),grady(1,2))
                     call addvecs(npts,gradz(1,3),gradz(1,1),gradz(1,2))

                           endif

      if((nfunctional.gt.0).and.(ipart.eq.0)) then

                     call addvecs(npts,hssxx(1,3),hssxx(1,1),hssxx(1,2))
                     call addvecs(npts,hssxy(1,3),hssxy(1,1),hssxy(1,2))
                     call addvecs(npts,hssxz(1,3),hssxz(1,1),hssxz(1,2))
                     call addvecs(npts,hssyy(1,3),hssyy(1,1),hssyy(1,2))
                     call addvecs(npts,hssyz(1,3),hssyz(1,1),hssyz(1,2))
                     call addvecs(npts,hsszz(1,3),hsszz(1,1),hsszz(1,2))

                                              endif

      do 1034 i=1,npts

      call xcfunctional(nfunctional,ixcfit,igrad,
     &                  dnsty(i,1),dnsty(i,2),dnsty(i,3),
     &                  gradx(i,1),grady(i,1),gradz(i,1),
     &                  gradx(i,2),grady(i,2),gradz(i,2),
     &                  gradx(i,3),grady(i,3),gradz(i,3),
     &                  hssxx(i,1),hssxy(i,1),hssxz(i,1),
     &                  hssyy(i,1),hssyz(i,1),hsszz(i,1),
     &                  hssxx(i,2),hssxy(i,2),hssxz(i,2),
     &                  hssyy(i,2),hssyz(i,2),hsszz(i,2),
     &                  hssxx(i,3),hssxy(i,3),hssxz(i,3),
     &                  hssyy(i,3),hssyz(i,3),hsszz(i,3),
     &                  xcnrg(i),xcpot(i,1),xcpot(i,2))

      if((nfunctional.gt.0).and.(ipart.eq.1)) then

          call xc_partial(nfunctional,xcnrg(i),xcpot(i,1),xcpot(i,2),
     &                    dnsty(i,1),dnsty(i,2),dnsty(i,3),
     &                    gradx(i,1),grady(i,1),gradz(i,1),
     &                    gradx(i,2),grady(i,2),gradz(i,2),
     &                    gradx(i,3),grady(i,3),gradz(i,3),
     &                    hssyy(i,3),hssyz(i,3),hsszz(i,3))

                                              endif

 1034 continue

      if((ixcfit.eq.0).and.(igrad.eq.0)) goto 1041

casa
c     if(igrad.eq.0) then

c                 do 1035 i=1,ncenters
c                 wkvec(i)=zero
c1035             if((dcc(i,iatom)-radius).le.extent4(i)) wkvec(i)=two

c                 do 1036 i=1,npts
c                 call partition1(ncenters,natomtype,xpoint(ja+i-1),
c    &                            ypoint(ja+i-1),zpoint(ja+i-1),wkvec,
c    &                            wkvec(ncenters+1),coord,
c    &                            weight_adjustment)

c                 xcpot(i,1)=xcpot(i,1)*weight_adjustment
c                 xcpot(i,2)=xcpot(i,2)*weight_adjustment

c1036             xcnrg(i)=xcnrg(i)*weight_adjustment

c                 call mpi_bsend(npts,1,
c    &                           mpi_integer,0,66,mpi_comm_world,ierr)

c                 call mpi_bsend(xcpot(1,1),npts,
c    &                           mpi_real8,0,66,mpi_comm_world,ierr)
c                 call mpi_bsend(xcpot(1,2),npts,
c    &                           mpi_real8,0,66,mpi_comm_world,ierr)
c                 call mpi_bsend(xcnrg(1),npts,
c    &                           mpi_real8,0,66,mpi_comm_world,ierr)

c                    endif

c     if((igrad.eq.0).and.(myid.eq.0)) then

c             do 1037 i=0,last_slave

c             call mpi_recv(nloop,1,
c    &                      mpi_integer,i,66,mpi_comm_world,istat,ierr)

c             call mpi_recv(wkvec(0*nloop+1),nloop,
c    &                      mpi_real8,i,66,mpi_comm_world,istat,ierr)
c             call mpi_recv(wkvec(1*nloop+1),nloop,
c    &                      mpi_real8,i,66,mpi_comm_world,istat,ierr)
c             call mpi_recv(wkvec(2*nloop+1),nloop,
c    &                      mpi_real8,i,66,mpi_comm_world,istat,ierr)

c             call fastwr(12,wkvec(0*nloop+1),nloop)
c             call fastwr(12,wkvec(1*nloop+1),nloop)
c             call fastwr(12,wkvec(2*nloop+1),nloop)

c1037         continue

c                                      endif

c1038 continue
casa

      if(igrad.eq.1) then
                           iwkvec(1)=ja
                           iwkvec(2)=npts
                           iwkvec(3)=myid

                           call mpi_bsend(iwkvec,3,mpi_integer,
     &                                    0,88,mpi_comm_world,ierr)

                           call mpi_bsend(xcpot(1,1),npts,mpi_real8,
     &                                    0,88,mpi_comm_world,ierr)
                           call mpi_bsend(xcpot(1,2),npts,mpi_real8,
     &                                    0,88,mpi_comm_world,ierr)
                           call mpi_bsend(xcnrg(1),npts,mpi_real8,
     &                                    0,88,mpi_comm_world,ierr)

                           goto 1074
casa
c                    else
c                          jhits=0

c                          do 1039 i=1,ncenters
c                          if(isys(i).eq.1) then

c                                 jhits=jhits+1
c                                 iwkvec(jhits)=i
c                                 distance=dcc(i,iatom)-radius
c                                 if(distance.le.extent4(i)) then

c                                                      wkvec(jhits)=two

c                                                            else

c                                                      wkvec(jhits)=zero

c                                                            endif

c                                             endif
c1039                      continue

c                          do 1040 i=1,npts
c                          call partition2(jhits,iwkvec,natomtype,
c    &                                     xpoint(ja+i-1),
c    &                                     ypoint(ja+i-1),
c    &                                     zpoint(ja+i-1),wkvec,
c    &                                     wkvec(ncenters+1),coord,
c    &                                     weight_adjustment)

c                          xcpot(i,1)=xcpot(i,1)*weight_adjustment
c                          xcpot(i,2)=xcpot(i,2)*weight_adjustment

c1040                      xcnrg(i)=xcnrg(i)*weight_adjustment
casaa
                     endif

 1041 continue

      jstart=ja-1

casa
c     if(ixcfit.eq.1) then
c                           do 1042 k=1,npts
c                           call weights(natomtype(iatom),
c    &                                   xpoint(jstart+k),
c    &                                   ypoint(jstart+k),
c    &                                   zpoint(jstart+k),
c    &                                   coord(1,iatom),adjust)
c1042                       wpoint(jstart+k)=wpoint(jstart+k)*adjust
c                     endif
casa

      do 1043 k=1,npts
      xcpot(k,1)=xcpot(k,1)*wpoint(jstart+k)
      xcpot(k,2)=xcpot(k,2)*wpoint(jstart+k)
 1043 xcnrg(k)=xcnrg(k)*wpoint(jstart+k)

      if(ixcfit.eq.0) then

             if(ipart.eq.1) then

                        do 1044 k=1,npts
                        hssyy(k,3)=hssyy(k,3)*wpoint(jstart+k)*two
                        hssyz(k,3)=hssyz(k,3)*wpoint(jstart+k)*two
 1044                   hsszz(k,3)=hsszz(k,3)*wpoint(jstart+k)

                        do 1045 k=1,npts
                        hssyy(k,1)=hssyy(k,3)*gradx(k,1)
                        hssyy(k,2)=hssyy(k,3)*grady(k,1)
 1045                   hssyy(k,3)=hssyy(k,3)*gradz(k,1)

                        do 1046 k=1,npts
                        hssyz(k,1)=hssyz(k,3)*gradx(k,2)
                        hssyz(k,2)=hssyz(k,3)*grady(k,2)
 1046                   hssyz(k,3)=hssyz(k,3)*gradz(k,2)

                        do 1047 k=1,npts
                        gradx(k,1)=hsszz(k,3)*gradx(k,1)
                        grady(k,1)=hsszz(k,3)*grady(k,1)
                        gradz(k,1)=hsszz(k,3)*gradz(k,1)
                        gradx(k,2)=hsszz(k,3)*gradx(k,2)
                        grady(k,2)=hsszz(k,3)*grady(k,2)
 1047                   gradz(k,2)=hsszz(k,3)*gradz(k,2)

                            endif

             ka=00*npts
             kb=01*npts
             kc=02*npts
             kd=03*npts
             ke=04*npts
             kf=05*npts
             kg=06*npts
             kh=07*npts
             ki=08*npts
             kj=09*npts
             kk=10*npts
             kl=11*npts
             km=12*npts
             kn=13*npts
             ko=14*npts
             kp=15*npts
             kq=16*npts
             kr=17*npts
             la=18*npts
             lb=19*npts
             lc=20*npts
             ld=21*npts
             le=22*npts
             lf=23*npts
             lg=24*npts
             lh=25*npts
             li=26*npts
             lj=27*npts
             lk=28*npts
             ll=29*npts

             ihits=0

             do 1048 i=1,ncontractions
             if(extent3(i).gt.tolerance) then

                        ihits=ihits+1
                        iwkvec(ihits)=i

                        ii=i*(i-1)/2

                        do 1049 k=1,npts
                        wkvec(ka+k)=fv(k,i)*xcpot(k,1)
                        wkvec(kb+k)=fv(k,i)*xcpot(k,2)
 1049                   wkvec(kc+k)=fv(k,i)*xcnrg(k)

                        if(nfunctional.eq.1) then

                                       do 1050 k=1,npts
                                       wkvec(kg+k)=fv(k,i)*hssyy(k,1)
                                       wkvec(kh+k)=fv(k,i)*hssyy(k,2)
                                       wkvec(ki+k)=fv(k,i)*hssyy(k,3)
                                       wkvec(kj+k)=fv(k,i)*hssyz(k,1)
                                       wkvec(kk+k)=fv(k,i)*hssyz(k,2)
                                       wkvec(kl+k)=fv(k,i)*hssyz(k,3)
                                       wkvec(km+k)=fvx(k,i)*hssyy(k,1)
                                       wkvec(kn+k)=fvy(k,i)*hssyy(k,2)
                                       wkvec(ko+k)=fvz(k,i)*hssyy(k,3)
                                       wkvec(kp+k)=fvx(k,i)*hssyz(k,1)
                                       wkvec(kq+k)=fvy(k,i)*hssyz(k,2)
                                       wkvec(kr+k)=fvz(k,i)*hssyz(k,3)
                                       wkvec(la+k)=fv(k,i)*gradx(k,2)
                                       wkvec(lb+k)=fv(k,i)*grady(k,2)
                                       wkvec(lc+k)=fv(k,i)*gradz(k,2)
                                       wkvec(ld+k)=fv(k,i)*gradx(k,1)
                                       wkvec(le+k)=fv(k,i)*grady(k,1)
                                       wkvec(lf+k)=fv(k,i)*gradz(k,1)
                                       wkvec(lg+k)=fvx(k,i)*gradx(k,2)
                                       wkvec(lh+k)=fvy(k,i)*grady(k,2)
                                       wkvec(li+k)=fvz(k,i)*gradz(k,2)
                                       wkvec(lj+k)=fvx(k,i)*gradx(k,1)
                                       wkvec(lk+k)=fvy(k,i)*grady(k,1)
 1050                                  wkvec(ll+k)=fvz(k,i)*gradz(k,1)

                                       do 1051 k=1,npts
                                       wkvec(km+k)=wkvec(km+k)
     &                                            +wkvec(kn+k)
     &                                            +wkvec(ko+k)
     &                                            +wkvec(lg+k)
     &                                            +wkvec(lh+k)
     &                                            +wkvec(li+k)
 1051                                  wkvec(kp+k)=wkvec(kp+k)
     &                                            +wkvec(kq+k)
     &                                            +wkvec(kr+k)
     &                                            +wkvec(lj+k)
     &                                            +wkvec(lk+k)
     &                                            +wkvec(ll+k)

                                       do 1052 k=1,npts
                                       wkvec(ke+k)=wkvec(km+k)
                                       wkvec(kf+k)=wkvec(kp+k)
                                       wkvec(kg+k)=wkvec(kg+k)
     &                                            +wkvec(la+k)
                                       wkvec(kh+k)=wkvec(kh+k)
     &                                            +wkvec(lb+k)
                                       wkvec(ki+k)=wkvec(ki+k)
     &                                            +wkvec(lc+k)
                                       wkvec(kj+k)=wkvec(kj+k)
     &                                            +wkvec(ld+k)
                                       wkvec(kk+k)=wkvec(kk+k)
     &                                            +wkvec(le+k)
 1052                                  wkvec(kl+k)=wkvec(kl+k)
     &                                            +wkvec(lf+k)

                                       do 1053 k=1,npts
                                       wkvec(ka+k)=wkvec(ka+k)
     &                                            +wkvec(ke+k)
 1053                                  wkvec(kb+k)=wkvec(kb+k)
     &                                            +wkvec(kf+k)

                                             endif

                       do 1054 j=1,ihits
                       jj=iwkvec(j)

                       if(extent3(i)*extent3(jj).lt.tolerance) goto 1054

                       ij=ii+jj

                       do 1055 k=1,npts
 1055                  focke(ij)=focke(ij)+fv(k,jj)*wkvec(kc+k)
                       do 1056 k=1,npts
 1056                  focka(ij)=focka(ij)+fv(k,jj)*wkvec(ka+k)
                       do 1057 k=1,npts
 1057                  fockb(ij)=fockb(ij)+fv(k,jj)*wkvec(kb+k)

                       if(nfunctional.eq.1) then

                               do 1058 k=1,npts
 1058                          focka(ij)=focka(ij)+fvx(k,jj)*wkvec(kg+k)
     &                                            +fvy(k,jj)*wkvec(kh+k)
     &                                            +fvz(k,jj)*wkvec(ki+k)
                               do 1059 k=1,npts
 1059                          fockb(ij)=fockb(ij)+fvx(k,jj)*wkvec(kj+k)
     &                                            +fvy(k,jj)*wkvec(kk+k)
     &                                            +fvz(k,jj)*wkvec(kl+k)

                                            endif

 1054                  continue

                                         endif
 1048        continue

                      endif

casa
c     if(ixcfit.eq.1) then

c             call gridxc(iatom,ncenters,npts,npoints,nxcfunctions,nxcs,
c    &                    nxcspd,ixclfunc,ixccfunc,radius,coord,alphaxc,
c    &                    coefsxc,coefpxc,coefdxc,xpoint(ja),ypoint(ja),
c    &                    zpoint(ja),fvx,extent2,extent,dcc,
c    &                    wkvec(0*npoints+1),wkvec(1*npoints+1),
c    &                    wkvec(2*npoints+1),wkvec(3*npoints+1),
c    &                    wkvec(4*npoints+1))

c             dlargest=dmaxvec(wpoint(ja),npts)

c             if(iteration.eq.1) then

c                         call sxcbld(ilast,nauxfunctions,npoints,
c    &                                nxcfunctions,npts,iwkvec,dlargest,
c    &                                tolerance,auxoverlap,extent,fvx,
c    &                                fvy,wpoint(ja),wkvec)

c                                endif

c             kstart=(nxcfunctions*(nxcfunctions+1))/2

c             do 1060 k=1,npts
c             if(dnsty(k,3).gt.small) then

c                      xcpot(k,1)=(xcpot(k,1)-xcnrg(k))/dnsty(k,3)
c                      xcpot(k,2)=(xcpot(k,2)-xcnrg(k))/dnsty(k,3)

c                                     else

c                      xcpot(k,1)=zero
c                      xcpot(k,2)=zero

c                                     endif
c1060         continue

c             do 1061 k=1,npts
c1061         wkvec(kstart+k)=zero

c             do 1062 i=1,nxcfunctions
c             if(extent(i).gt.tolerance) then

c                      do 1063 k=1,npts
c1063                  wkvec(kstart+k)=wkvec(kstart+k)+tvecb(i)*fvx(k,i)

c                                        endif
c1062         continue

c             ka=1*npts+kstart
c             kb=2*npts+kstart

c             do 1064 k=1,npts
c1064         wkvec(ka+k)=wkvec(kstart+k)*xcpot(k,1)
c             do 1065 k=1,npts
c1065         wkvec(kb+k)=wkvec(kstart+k)*xcpot(k,2)

c             ihits=0

c             do 1066 i=1,ncontractions
c             if(extent3(i).gt.tolerance) then

c                      ihits=ihits+1
c                      iwkvec(ihits)=i

c                      ii=i*(i-1)/2

c                      kc=3*npts+kstart
c                      kd=4*npts+kstart

c                      do 1067 k=1,npts
c1067                  wkvec(kc+k)=fv(k,i)*wkvec(ka+k)
c                      do 1068 k=1,npts
c1068                  wkvec(kd+k)=fv(k,i)*wkvec(kb+k)

c                      do 1069 j=1,ihits
c                      jj=iwkvec(j)
c                      if(extent3(i)*extent3(jj).lt.tolerance) goto 1069
c                      ij=ii+jj
c                      do 1070 k=1,npts
c1070                  focka(ij)=focka(ij)+fv(k,jj)*wkvec(kc+k)
c                      do 1071 k=1,npts
c1071                  fockb(ij)=fockb(ij)+fv(k,jj)*wkvec(kd+k)
c1069                  continue

c                                         endif
c1066 continue
c             
c                     endif

c     do 1072 i=1,nxcfunctions
c     if(extent(i)*dlargest.gt.tolerance) then

c                            do 1073 k=1,npts
c1073                        tvece(i)=tvece(i)+fvx(k,i)*xcnrg(k)

c                                         endif
c1072 continue
casa

 1074 continue

      if((myid.eq.0).and.(igrad.eq.1)) then

              do 1075 i=1,last_slave+1

              call mpi_recv(iwkvec,3,mpi_integer,
     &                      mpi_any_source,88,mpi_comm_world,istat,ierr)

              ja=iwkvec(1)
              npts=iwkvec(2)
              islave=iwkvec(3)

              call mpi_recv(xpoint(ja),npts,mpi_real8,
     &                      islave,88,mpi_comm_world,istat,ierr)
              call mpi_recv(ypoint(ja),npts,mpi_real8,
     &                      islave,88,mpi_comm_world,istat,ierr)
              call mpi_recv(zpoint(ja),npts,mpi_real8,
     &                      islave,88,mpi_comm_world,istat,ierr)

 1075         continue

                                       endif

 1005 continue

      if((myid.eq.0).and.(igrad.eq.1)) then

                                call fastwr(12,xpoint,nptsatom(iatom))
                                call fastwr(12,ypoint,nptsatom(iatom))
                                call fastwr(12,zpoint,nptsatom(iatom))

                                       endif

 1004 continue

casa
c     if(ixcfit.eq.1) then

c             if(iteration.eq.1) then

c                     na=nxcfunctions*(nxcfunctions+1)/2
c
c                     call mpi_reduce(wkvec(1),wkvec(1+na),na,mpi_real8,
c    &                                mpi_sum,0,mpi_comm_world,ierr)
c             
c                     if(myid.eq.0) call overlpcopy(nauxfunctions,
c    &                                              nxcfunctions,
c    &                                              wkvec(na+1),
c    &                                              auxoverlap)

c                                endif

c             call mpi_reduce(tvece,wkvec,nxcfunctions,
c    &                        mpi_real8,mpi_sum,0,mpi_comm_world,ierr)
c
c             if(myid.eq.0) call stringcopy(nxcfunctions,tvece,wkvec)
c
c             call mpi_bcast(tvece,nxcfunctions,
c    &                       mpi_real8,0,mpi_comm_world,ierr)

c                     endif
casa

      nfunctional=nsavefunctional

      return
      end
