      subroutine gradmm(nreal,ninteger,nfockdim,ncenters,nconts,ncontp,
     &                  ncontd,ncds,ncdspd,nshels,nshelp,nsheld,mtloca,
     &                  iprimvc,nprimvc,icentvc,icdlfunc,icdcfunc,
     &                  iwkvec,coord,charge,alpha,coeff,cdfitc,alphacd,
     &                  coefscd,coefpcd,coefdcd,gradient,overlap,dmata,
     &                  dmatb,density,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 jsend(30)

      dimension istat(mpi_status_size)

      dimension g(8)

      dimension nshels(*),nshelp(*),nsheld(*),mtloca(*),iprimvc(*),
     &          nprimvc(*),icentvc(*),icdlfunc(*),icdcfunc(*),iwkvec(*)

      dimension coord(3,*),charge(*),alpha(*),coeff(*),cdfitc(*),
     &          alphacd(*),coefscd(*),coefpcd(*),coefdcd(*),
     &          gradient(3,*),overlap(nfockdim,*),dmata(*),dmatb(*),
     &          density(*),wkvec(*)

      data zero,half,one,three/0.0d0,0.5d0,1.0d0,3.0d0/

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

      call igftable

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

      call mpi_bcast(charge,ncenters,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(dmata,nfockdim,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(dmatb,nfockdim,mpi_real8,0,mpi_comm_world,ierr)

      sqrt3=sqrt(three)

      do 1001 i=1,nfockdim
 1001 density(i)=dmata(i)+dmatb(i)

      do 1002 i=1,ncontd
      do 1002 j=1,nconts

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshels(j))
      i4=i2+i3

      density(i4)=density(i4)/sqrt3
      density(i4+i1*3+3)=density(i4+i1*3+3)/sqrt3
      density(i4+i1*5+10)=density(i4+i1*5+10)/sqrt3

 1002 continue

      do 1003 i=1,ncontd
      do 1003 j=1,ncontp

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nshelp(j))
      i4=i2+i3

      density(i4)=density(i4)/sqrt3
      density(i4+1)=density(i4+1)/sqrt3
      density(i4+2)=density(i4+2)/sqrt3

      density(i4+i1*3+3)=density(i4+i1*3+3)/sqrt3
      density(i4+i1*3+4)=density(i4+i1*3+4)/sqrt3
      density(i4+i1*3+5)=density(i4+i1*3+5)/sqrt3

      density(i4+i1*5+10)=density(i4+i1*5+10)/sqrt3
      density(i4+i1*5+11)=density(i4+i1*5+11)/sqrt3
      density(i4+i1*5+12)=density(i4+i1*5+12)/sqrt3

 1003 continue

      do 1004 i=1,ncontd
      do 1004 j=1,i

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nsheld(j))
      i4=i2+i3

      density(i4)=density(i4)/three
      if(i.ne.j) density(i4+1)=density(i4+1)/sqrt3
      if(i.ne.j) density(i4+2)=density(i4+2)/sqrt3
      if(i.ne.j) density(i4+3)=density(i4+3)/three
      if(i.ne.j) density(i4+4)=density(i4+4)/sqrt3
      if(i.ne.j) density(i4+5)=density(i4+5)/three

      density(i4+i1)=density(i4+i1)/sqrt3
      if(i.ne.j) density(i4+i1+3)=density(i4+i1+3)/sqrt3
      if(i.ne.j) density(i4+i1+5)=density(i4+i1+5)/sqrt3

      density(i4+i1*2+1)=density(i4+i1*2+1)/sqrt3
      if(i.ne.j) density(i4+i1*2+4)=density(i4+i1*2+4)/sqrt3
      if(i.ne.j) density(i4+i1*2+6)=density(i4+i1*2+6)/sqrt3

      density(i4+i1*3+3)=density(i4+i1*3+3)/three
      density(i4+i1*3+4)=density(i4+i1*3+4)/sqrt3
      density(i4+i1*3+5)=density(i4+i1*3+5)/sqrt3
      density(i4+i1*3+6)=density(i4+i1*3+6)/three
      if(i.ne.j) density(i4+i1*3+7)=density(i4+i1*3+7)/sqrt3
      if(i.ne.j) density(i4+i1*3+8)=density(i4+i1*3+8)/three

      density(i4+i1*4+6)=density(i4+i1*4+6)/sqrt3
      density(i4+i1*4+9)=density(i4+i1*4+9)/sqrt3
      if(i.ne.j) density(i4+i1*4+11)=density(i4+i1*4+11)/sqrt3

      density(i4+i1*5+10)=density(i4+i1*5+10)/three
      density(i4+i1*5+11)=density(i4+i1*5+11)/sqrt3
      density(i4+i1*5+12)=density(i4+i1*5+12)/sqrt3
      density(i4+i1*5+13)=density(i4+i1*5+13)/three
      density(i4+i1*5+14)=density(i4+i1*5+14)/sqrt3
      density(i4+i1*5+15)=density(i4+i1*5+15)/three

 1004 continue

      do 1005 i=1,ncenters
      do 1005 j=1,3
 1005 gradient(j,i)=zero

      do 1006 j=1,3
      do 1006 i=1,nfockdim
 1006 overlap(i,j)=zero

      if(myid.eq.0) then

              do 1007 i=1,ncenters
              do 1007 j=1,i

              if(i.eq.j) goto 1007

              factor=charge(i)*charge(j)/
     &               ((((coord(1,i)-coord(1,j))**2)+
     &                 ((coord(2,i)-coord(2,j))**2)+
     &                 ((coord(3,i)-coord(3,j))**2))**(one+half))

              gradient(1,i)=gradient(1,i)+factor*(coord(1,j)-coord(1,i))
              gradient(2,i)=gradient(2,i)+factor*(coord(2,j)-coord(2,i))
              gradient(3,i)=gradient(3,i)+factor*(coord(3,j)-coord(3,i))

              gradient(1,j)=gradient(1,j)+factor*(coord(1,i)-coord(1,j))
              gradient(2,j)=gradient(2,j)+factor*(coord(2,i)-coord(2,j))
              gradient(3,j)=gradient(3,j)+factor*(coord(3,i)-coord(3,j))

 1007         continue

                      endif

      ncposs=nconts*(nconts+1)/2
      ncpopp=ncontp*(ncontp+1)/2
      ncpodd=ncontd*(ncontd+1)/2

      ncpops=ncontp*nconts
      ncpods=ncontd*nconts
      ncpodp=ncontd*ncontp

      npposs=0
      nppops=0
      nppopp=0
      nppods=0
      nppodp=0
      nppodd=0

      do 1008 i=1,nconts
      do 1008 j=1,i
 1008 npposs=npposs+nprimvc(nshels(i))*nprimvc(nshels(j))

      do 1009 i=1,ncontp
      do 1009 j=1,nconts
 1009 nppops=nppops+nprimvc(nshelp(i))*nprimvc(nshels(j))

      do 1010 i=1,ncontp
      do 1010 j=1,i
 1010 nppopp=nppopp+nprimvc(nshelp(i))*nprimvc(nshelp(j))

      do 1011 i=1,ncontd
      do 1011 j=1,nconts
 1011 nppods=nppods+nprimvc(nsheld(i))*nprimvc(nshels(j))

      do 1012 i=1,ncontd
      do 1012 j=1,ncontp
 1012 nppodp=nppodp+nprimvc(nsheld(i))*nprimvc(nshelp(j))

      do 1013 i=1,ncontd
      do 1013 j=1,i
 1013 nppodd=nppodd+nprimvc(nsheld(i))*nprimvc(nsheld(j))

      ja=1
      jb=ja+npposs
      jc=jb+ncposs

      ia=jc+ncposs*2
      ib=ia+npposs

      if(ib.gt.ninteger) stop 'gradmm ss integer overflow'

      je=1
      jf=je+npposs
      jg=jf+npposs
      jh=jg+npposs
      ji=jh+npposs
      jj=ji+npposs
      jk=jj+npposs
      jl=jk+npposs
      jm=jl+npposs
      jn=jm+npposs
      jo=jn+npposs
      jp=jo+npposs
      jq=jp+npposs
      jr=jq+npposs
      js=jr+npposs
      jt=js+npposs
      ju=jt+npposs
      jv=ju+npposs
      jw=jv+npposs
      jx=jw+ncposs
      jy=jx+ncposs
      jz=jy+ncposs
      j0=jz+npposs*13

      if(j0.gt.nreal) stop 'gradmm ss real overflow'

      call setupgss(ncposs,n,nconts,nshels,mtloca,iprimvc,nprimvc,
     &              icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),coord,
     &              alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),wkvec(jh),
     &              wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),wkvec(jm),
     &              wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),wkvec(jr),
     &              wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jw),
     &              wkvec(jx),wkvec(jy),wkvec(jz))

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                          call stringcopy(n,wkvec(15*n+1),wkvec(jt))
                          call stringcopy(n,wkvec(16*n+1),wkvec(ju))
                          call stringcopy(n,wkvec(17*n+1),wkvec(jv))
                    endif

      call stringcopy(ncposs,wkvec(18*n+0*ncposs+1),wkvec(jw))
      call stringcopy(ncposs,wkvec(18*n+1*ncposs+1),wkvec(jx))
      call stringcopy(ncposs,wkvec(18*n+2*ncposs+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncposs
      jy=jx+ncposs
      jz=jy+ncposs

      if(myid.eq.0) then
                          isent=0

                          do 1014 islave=0,numprocs-1

                          i1=n-isent
                          i2=numprocs-islave

                          nloop=i1/i2
                          if(mod(i1,i2).ne.0) nloop=nloop+1

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

                          jsend(01)=ja
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js
                          jsend(17)=jt
                          jsend(18)=ju
                          jsend(19)=jv

                          call mpi_bsend(jsend,19,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

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

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jh),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ji),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jt),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ju),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jv),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          ja=ja+nloop
 
                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop
                          jt=jt+nloop
                          ju=ju+nloop
                          jv=jv+nloop
 
 1014                     isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,19,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ja=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)
      jt=jsend(17)
      ju=jsend(18)
      jv=jsend(19)

      call mpi_recv(iwkvec(ja),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jh),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ji),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jt),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ju),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jv),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop*3
      kd=kc+ncposs*3
      ke=kd+ncposs*3

      if(ke.gt.nreal) stop 'gradmm ss real overflow'

      call gradovss(nfockdim,ncposs,nloop,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),density,overlap,gradient,wkvec(je),
     &              wkvec(jh),wkvec(jm),wkvec(jn),wkvec(jp),wkvec(jo),
     &              wkvec(jt),wkvec(ju),wkvec(jv),wkvec(ka),wkvec(kb),
     &              wkvec(kc),wkvec(kd))

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop*3
      ki=kh+ncposs
      kj=ki+ncposs*3
      kk=kj+ncposs*3
      kl=kk+ncposs*3

      if(kl.gt.nreal) stop 'gradmm ss real overflow'

      call gradnass(ncposs,nloop,ncenters,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),g,coord,charge,density,gradient,
     &              wkvec(je),wkvec(jh),wkvec(ji),wkvec(jq),wkvec(jr),
     &              wkvec(js),wkvec(jp),wkvec(jt),wkvec(ju),wkvec(jv),
     &              wkvec(jw),wkvec(jx),wkvec(jy),wkvec(ka),wkvec(kb),
     &              wkvec(kc),wkvec(kd),wkvec(ke),wkvec(kf),wkvec(kg),
     &              wkvec(kh),wkvec(ki),wkvec(kj),wkvec(kk))

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop*3
      ks=kr+ncposs
      kt=ks+ncposs*3
      ku=kt+ncposs*3
      kv=ku+ncposs*3

      if(kv.gt.nreal) stop 'gradmm ss real overflow'

      call gradcdss(ncposs,nloop,ncds,icdlfunc,icdcfunc,iwkvec(ja),
     &              iwkvec(jb),iwkvec(jc),iwkvec(ia),g,coord,density,
     &              gradient,cdfitc,alphacd,coefscd,wkvec(je),
     &              wkvec(jh),wkvec(ji),wkvec(jq),wkvec(jr),wkvec(js),
     &              wkvec(jp),wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jw),
     &              wkvec(jx),wkvec(jy),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))

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop
      kv=ku+nloop
      kw=kv+nloop
      kx=kw+nloop
      ky=kx+nloop*3
      kz=ky+nloop*3
      la=kz+nloop*3
      lb=la+nloop*3
      lc=lb+nloop*3
      ld=lc+nloop*6
      le=ld+nloop*9
      lf=le+nloop*9
      lg=lf+nloop*18
      lh=lg+nloop
      li=lh+nloop*3
      lj=li+nloop*3
      lk=lj+ncposs
      ll=lk+ncposs*3
      lm=ll+ncposs*6
      ln=lm+ncposs*3
      lo=ln+ncposs*9
      lp=lo+ncposs*18
      lq=lp+ncposs*3
      lr=lq+ncposs*9
      ls=lr+ncposs*18
      lt=ls+ncposs*3

      if(lt.gt.nreal) stop 'gradmm ss real overflow'

      call gradxxss(ncposs,nloop,ncds,ncdspd,icdlfunc,icdcfunc,
     &              iwkvec(ja),iwkvec(jb),iwkvec(jc),iwkvec(ia),g,coord,
     &              density,gradient,cdfitc,alphacd,coefscd,coefpcd,
     &              coefdcd,wkvec(je),wkvec(jh),wkvec(ji),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jp),wkvec(jt),wkvec(ju),
     &              wkvec(jv),wkvec(jw),wkvec(jx),wkvec(jy),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),wkvec(ln),wkvec(lo),
     &              wkvec(lp),wkvec(lq),wkvec(lr),wkvec(ls))

      if(ncontp.eq.0) goto 2001

      ja=1
      jb=ja+nppops
      jc=jb+ncpops*3

      ia=jc+ncpops*2
      ib=ia+nppops

      if(ib.gt.ninteger) stop 'gradmm ps integer overflow'

      je=1
      jf=je+nppops
      jg=jf+nppops
      jh=jg+nppops
      ji=jh+nppops
      jj=ji+nppops
      jk=jj+nppops
      jl=jk+nppops
      jm=jl+nppops
      jn=jm+nppops
      jo=jn+nppops
      jp=jo+nppops
      jq=jp+nppops
      jr=jq+nppops
      js=jr+nppops
      jt=js+nppops
      ju=jt+nppops
      jv=ju+nppops
      jw=jv+nppops
      jx=jw+ncpops
      jy=jx+ncpops
      jz=jy+ncpops
      j0=jz+nppops*13

      if(j0.gt.nreal) stop 'gradmm ps real overflow'

      call setupgps(ncpops,n,ncontp,nconts,nshelp,nshels,mtloca,iprimvc,
     &              nprimvc,icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),
     &              coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &              wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &              wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &              wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jz))

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                          call stringcopy(n,wkvec(15*n+1),wkvec(jt))
                          call stringcopy(n,wkvec(16*n+1),wkvec(ju))
                          call stringcopy(n,wkvec(17*n+1),wkvec(jv))
                    endif

      call stringcopy(ncpops,wkvec(18*n+0*ncpops+1),wkvec(jw))
      call stringcopy(ncpops,wkvec(18*n+1*ncpops+1),wkvec(jx))
      call stringcopy(ncpops,wkvec(18*n+2*ncpops+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpops
      jy=jx+ncpops
      jz=jy+ncpops

      if(myid.eq.0) then
                          isent=0

                          do 1015 islave=0,numprocs-1

                          i1=n-isent
                          i2=numprocs-islave

                          nloop=i1/i2
                          if(mod(i1,i2).ne.0) nloop=nloop+1

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

                          jsend(01)=ja
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js
                          jsend(17)=jt
                          jsend(18)=ju
                          jsend(19)=jv

                          call mpi_bsend(jsend,19,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

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

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jh),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ji),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jt),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ju),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jv),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          ja=ja+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop
                          jt=jt+nloop
                          ju=ju+nloop
                          jv=jv+nloop

1015                      isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,19,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ja=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)
      jt=jsend(17)
      ju=jsend(18)
      jv=jsend(19)

      call mpi_recv(iwkvec(ja),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jh),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ji),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jt),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ju),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jv),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop*3
      kf=ke+nloop*6
      kg=kf+nloop
      kh=kg+ncpops
      ki=kh+ncpops*6
      kj=ki+ncpops*9

      if(kj.gt.nreal) stop 'gradmm ps real overflow'

      call gradovps(nfockdim,ncpops,nloop,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),density,overlap,gradient,wkvec(je),
     &              wkvec(jh),wkvec(jm),wkvec(jn),wkvec(jp),wkvec(jo),
     &              wkvec(jt),wkvec(ju),wkvec(jv),wkvec(ka),wkvec(kb),
     &              wkvec(kc),wkvec(kd),wkvec(ke),wkvec(kf),wkvec(kg),
     &              wkvec(kh),wkvec(ki))

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop*3
      kk=kj+nloop*3
      kl=kk+nloop*6
      km=kl+nloop
      kn=km+ncpops
      ko=kn+ncpops*3
      kp=ko+ncpops*6
      kq=kp+ncpops*9
      kr=kq+ncpops*9

      if(kr.gt.nreal) stop 'gradmm ps real overflow'

      call gradnaps(ncpops,nloop,ncenters,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),g,coord,charge,density,gradient,
     &              wkvec(je),wkvec(jh),wkvec(ji),wkvec(jq),wkvec(jr),
     &              wkvec(js),wkvec(jp),wkvec(jt),wkvec(ju),wkvec(jv),
     &              wkvec(jw),wkvec(jx),wkvec(jy),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))

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop*3
      ku=kt+nloop*3
      kv=ku+nloop*6
      kw=kv+nloop
      kx=kw+ncpops
      ky=kx+ncpops*3
      kz=ky+ncpops*6
      k0=kz+ncpops*9
      k1=k0+ncpops*9

      if(k1.gt.nreal) stop 'gradmm ps real overflow'

      call gradcdps(ncpops,nloop,ncds,icdlfunc,icdcfunc,iwkvec(ja),
     &              iwkvec(jb),iwkvec(jc),iwkvec(ia),g,coord,density,
     &              gradient,cdfitc,alphacd,coefscd,wkvec(je),wkvec(jh),
     &              wkvec(ji),wkvec(jq),wkvec(jr),wkvec(js),wkvec(jp),
     &              wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jw),wkvec(jx),
     &              wkvec(jy),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(k0))

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop
      kv=ku+nloop
      kw=kv+nloop
      kx=kw+nloop
      ky=kx+nloop
      kz=ky+nloop*3
      la=kz+nloop*3
      lb=la+nloop*3
      lc=lb+nloop*3
      ld=lc+nloop*6
      le=ld+nloop*6
      lf=le+nloop*6
      lg=lf+nloop*3
      lh=lg+nloop*3
      li=lh+nloop*6
      lj=li+nloop*9
      lk=lj+nloop*9
      ll=lk+nloop*18
      lm=ll+nloop*18
      ln=lm+nloop*18
      lo=ln+nloop*36
      lp=lo+nloop
      lq=lp+nloop*6
      lr=lq+nloop*9
      ls=lr+ncpops
      lt=ls+ncpops*3
      lu=lt+ncpops*6
      lv=lu+ncpops*3
      lw=lv+ncpops*9
      lx=lw+ncpops*18
      ly=lx+ncpops*6
      lz=ly+ncpops*18
      ma=lz+ncpops*36
      mb=ma+ncpops*9
      mc=mb+ncpops*27
      md=mc+ncpops*54
      me=md+ncpops*9

      if(me.gt.nreal) stop 'gradmm ps real overflow'

      call gradxxps(ncpops,nloop,ncds,ncdspd,icdlfunc,icdcfunc,
     &              iwkvec(ja),iwkvec(jb),iwkvec(jc),iwkvec(ia),g,coord,
     &              density,gradient,cdfitc,alphacd,coefscd,coefpcd,
     &              coefdcd,wkvec(je),wkvec(jh),wkvec(ji),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jp),wkvec(jt),wkvec(ju),
     &              wkvec(jv),wkvec(jw),wkvec(jx),wkvec(jy),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),wkvec(ln),wkvec(lo),
     &              wkvec(lp),wkvec(lq),wkvec(lr),wkvec(ls),wkvec(lt),
     &              wkvec(lu),wkvec(lv),wkvec(lw),wkvec(lx),wkvec(ly),
     &              wkvec(lz),wkvec(ma),wkvec(mb),wkvec(mc),wkvec(md))

      ja=1
      jb=ja+ncpopp
      jc=jb+nppopp
      jd=jc+ncpopp*9

      ia=jd+ncpopp*2
      ib=ia+nppopp

      if(ib.gt.ninteger) stop 'gradmm pp integer overflow'

      je=1
      jf=je+nppopp
      jg=jf+nppopp
      jh=jg+nppopp
      ji=jh+nppopp
      jj=ji+nppopp
      jk=jj+nppopp
      jl=jk+nppopp
      jm=jl+nppopp
      jn=jm+nppopp
      jo=jn+nppopp
      jp=jo+nppopp
      jq=jp+nppopp
      jr=jq+nppopp
      js=jr+nppopp
      jt=js+nppopp
      ju=jt+nppopp
      jv=ju+nppopp
      jw=jv+nppopp
      jx=jw+nppopp
      jy=jx+nppopp
      jz=jy+nppopp
      j0=jz+ncpopp
      j1=j0+ncpopp
      j2=j1+ncpopp
      j3=j2+nppopp*16

      if(j3.gt.nreal) stop 'gradmm pp real overflow'

      call setupgpp(ncpopp,n,ncontp,nshelp,mtloca,iprimvc,nprimvc,
     &              icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),iwkvec(jd),
     &              coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &              wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &              wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &              wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jz),wkvec(j0),
     &              wkvec(j1),wkvec(j2))


      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                          call stringcopy(n,wkvec(15*n+1),wkvec(jt))
                          call stringcopy(n,wkvec(16*n+1),wkvec(ju))
                          call stringcopy(n,wkvec(17*n+1),wkvec(jv))
                          call stringcopy(n,wkvec(18*n+1),wkvec(jw))
                          call stringcopy(n,wkvec(19*n+1),wkvec(jx))
                          call stringcopy(n,wkvec(20*n+1),wkvec(jy))
                    endif

      call stringcopy(ncpopp,wkvec(21*n+0*ncpopp+1),wkvec(jz))
      call stringcopy(ncpopp,wkvec(21*n+1*ncpopp+1),wkvec(j0))
      call stringcopy(ncpopp,wkvec(21*n+2*ncpopp+1),wkvec(j1))
 
      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+n
      jy=jx+n
      jz=jy+n
      j0=jz+ncpopp
      j1=j0+ncpopp
      j2=j1+ncpopp

      if(myid.eq.0) then
                          isent=0

                          do 1016 islave=0,numprocs-1

                          i1=n-isent
                          i2=numprocs-islave

                          nloop=i1/i2
                          if(mod(i1,i2).ne.0) nloop=nloop+1

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

                          jsend(01)=jb
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js
                          jsend(17)=jt
                          jsend(18)=ju
                          jsend(19)=jv
                          jsend(20)=jw
                          jsend(21)=jx
                          jsend(22)=jy

                          call mpi_bsend(jsend,22,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          call mpi_bsend(iwkvec(jb),nloop,
     &                                   mpi_integer,islave,islave,
     &                                   mpi_comm_world,ierr)

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jh),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ji),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jt),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ju),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jv),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jw),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jx),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jy),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          jb=jb+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop
                          jt=jt+nloop
                          ju=ju+nloop
                          jv=jv+nloop
                          jw=jw+nloop
                          jx=jx+nloop
                          jy=jy+nloop

1016                      isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,22,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      jb=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)
      jt=jsend(17)
      ju=jsend(18)
      jv=jsend(19)
      jw=jsend(20)
      jx=jsend(21)
      jy=jsend(22)

      call mpi_recv(iwkvec(jb),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jh),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ji),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jt),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ju),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jv),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jw),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jx),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jy),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop*3
      kf=ke+nloop*3
      kg=kf+nloop*6
      kh=kg+nloop*18
      ki=kh+nloop
      kj=ki+nloop*3
      kk=kj+ncpopp*3
      kl=kk+ncpopp*18
      km=kl+ncpopp*27

      if(km.gt.nreal) stop 'gradmm pp real overflow'

      call gradovpp(nfockdim,ncpopp,nloop,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),iwkvec(jd),density,overlap,gradient,
     &              wkvec(je),wkvec(jh),wkvec(jm),wkvec(jn),wkvec(jp),
     &              wkvec(jo),wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jw),
     &              wkvec(jx),wkvec(jy),wkvec(ka),wkvec(kb),wkvec(kc),
     &              wkvec(kd),wkvec(ke),wkvec(kf),wkvec(kg),wkvec(kh),
     &              wkvec(ki),wkvec(kj),wkvec(kk),wkvec(kl))

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop*3
      kl=kk+nloop*3
      km=kl+nloop*3
      kn=km+nloop*6
      ko=kn+nloop*6
      kp=ko+nloop*10
      kq=kp+nloop
      kr=kq+nloop*3
      ks=kr+ncpopp
      kt=ks+ncpopp*3
      ku=kt+ncpopp*3
      kv=ku+ncpopp*6
      kw=kv+ncpopp*9
      kx=kw+ncpopp*10
      ky=kx+ncpopp*18
      kz=ky+ncpopp*18
      k0=kz+ncpopp*27

      if(k0.gt.nreal) stop 'gradmm pp real overflow'

      call gradnapp(ncpopp,nloop,ncenters,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),iwkvec(jd),g,coord,charge,density,
     &              gradient,wkvec(je),wkvec(jh),wkvec(ji),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jp),wkvec(jt),wkvec(ju),
     &              wkvec(jv),wkvec(jz),wkvec(j0),wkvec(j1),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))

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop 
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop*3
      kv=ku+nloop*3
      kw=kv+nloop*3
      kx=kw+nloop*6
      ky=kx+nloop*6
      kz=ky+nloop*10
      la=kz+nloop
      lb=la+nloop*3
      lc=lb+ncpopp
      ld=lc+ncpopp*3
      le=ld+ncpopp*3
      lf=le+ncpopp*6
      lg=lf+ncpopp*9
      lh=lg+ncpopp*10
      li=lh+ncpopp*18
      lj=li+ncpopp*18
      lk=lj+ncpopp*27

      if(lk.gt.nreal) stop 'gradmm pp real overflow'

      call gradcdpp(ncpopp,nloop,ncds,icdlfunc,icdcfunc,iwkvec(ja),
     &              iwkvec(jb),iwkvec(jc),iwkvec(jd),iwkvec(ia),g,coord,
     &              density,gradient,cdfitc,alphacd,coefscd,wkvec(je),
     &              wkvec(jh),wkvec(ji),wkvec(jq),wkvec(jr),wkvec(js),
     &              wkvec(jp),wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jz),
     &              wkvec(j0),wkvec(j1),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))

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop
      kv=ku+nloop
      kw=kv+nloop
      kx=kw+nloop
      ky=kx+nloop
      kz=ky+nloop
      la=kz+nloop*3
      lb=la+nloop*3
      lc=lb+nloop*3
      ld=lc+nloop*3
      le=ld+nloop*3
      lf=le+nloop*6
      lg=lf+nloop*6
      lh=lg+nloop*6
      li=lh+nloop*6
      lj=li+nloop*10
      lk=lj+nloop*10
      ll=lk+nloop*10
      lm=ll+nloop*3
      ln=lm+nloop*3
      lo=ln+nloop*6
      lp=lo+nloop*9
      lq=lp+nloop*9
      lr=lq+nloop*18
      ls=lr+nloop*18
      lt=ls+nloop*18
      lu=lt+nloop*36
      lv=lu+nloop*30
      lw=lv+nloop*30
      lx=lw+nloop*60
      ly=lx+nloop
      lz=ly+nloop*10
      ma=lz+nloop*18
      mb=ma+ncpopp
      mc=mb+ncpopp*3
      md=mc+ncpopp*6
      me=md+ncpopp*3
      mf=me+ncpopp*9
      mg=mf+ncpopp*18
      mh=mg+ncpopp*3
      mi=mh+ncpopp*9
      mj=mi+ncpopp*18
      mk=mj+ncpopp*6
      ml=mk+ncpopp*18
      mm=ml+ncpopp*36
      mn=mm+ncpopp*9
      mo=mn+ncpopp*27
      mp=mo+ncpopp*54
      mq=mp+ncpopp*10
      mr=mq+ncpopp*30
      ms=mr+ncpopp*60
      mt=ms+ncpopp*18
      mu=mt+ncpopp*54
      mv=mu+ncpopp*108
      mw=mv+ncpopp*18
      mx=mw+ncpopp*54
      my=mx+ncpopp*108
      mz=my+ncpopp*27

      if(mz.gt.nreal) stop 'gradmm pp real overflow'

      call gradxxpp(ncpopp,nloop,ncds,ncdspd,icdlfunc,icdcfunc,
     &              iwkvec(ja),iwkvec(jb),iwkvec(jc),iwkvec(jd),
     &              iwkvec(ia),g,coord,density,gradient,cdfitc,alphacd,
     &              coefscd,coefpcd,coefdcd,wkvec(je),wkvec(jh),
     &              wkvec(ji),wkvec(jq),wkvec(jr),wkvec(js),wkvec(jp),
     &              wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jz),wkvec(j0),
     &              wkvec(j1),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),
     &              wkvec(ln),wkvec(lo),wkvec(lp),wkvec(lq),wkvec(lr),
     &              wkvec(ls),wkvec(lt),wkvec(lu),wkvec(lv),wkvec(lw),
     &              wkvec(lx),wkvec(ly),wkvec(lz),wkvec(ma),wkvec(mb),
     &              wkvec(mc),wkvec(md),wkvec(me),wkvec(mf),wkvec(mg),
     &              wkvec(mh),wkvec(mi),wkvec(mj),wkvec(mk),wkvec(ml),
     &              wkvec(mm),wkvec(mn),wkvec(mo),wkvec(mp),wkvec(mq),
     &              wkvec(mr),wkvec(ms),wkvec(mt),wkvec(mu),wkvec(mv),
     &              wkvec(mw),wkvec(mx),wkvec(my))

      if(ncontd.eq.0) goto 2001

      ja=1
      jb=ja+nppods
      jc=jb+ncpods*6

      ia=jc+ncpods*2
      ib=ia+nppods

      if(ib.gt.ninteger) stop 'gradmm ds integer overflow'

      je=1
      jf=je+nppods
      jg=jf+nppods
      jh=jg+nppods
      ji=jh+nppods
      jj=ji+nppods
      jk=jj+nppods
      jl=jk+nppods
      jm=jl+nppods
      jn=jm+nppods
      jo=jn+nppods
      jp=jo+nppods
      jq=jp+nppods
      jr=jq+nppods
      js=jr+nppods
      jt=js+nppods
      ju=jt+nppods
      jv=ju+nppods
      jw=jv+nppods
      jx=jw+ncpods
      jy=jx+ncpods
      jz=jy+ncpods
      j1=jz+nppods*13

      if(j1.gt.nreal) stop 'gradmm ds real overflow'

      call setupgds(ncpods,n,ncontd,nconts,nsheld,nshels,mtloca,iprimvc,
     &              nprimvc,icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),
     &              coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &              wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &              wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &              wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jz))

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                          call stringcopy(n,wkvec(15*n+1),wkvec(jt))
                          call stringcopy(n,wkvec(16*n+1),wkvec(ju))
                          call stringcopy(n,wkvec(17*n+1),wkvec(jv))
                    endif
 
      call stringcopy(ncpods,wkvec(18*n+0*ncpods+1),wkvec(jw))
      call stringcopy(ncpods,wkvec(18*n+1*ncpods+1),wkvec(jx))
      call stringcopy(ncpods,wkvec(18*n+2*ncpods+1),wkvec(jy))
 
      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpods
      jy=jx+ncpods
      jz=jy+ncpods

      if(myid.eq.0) then
                          isent=0

                          do 1017 islave=0,numprocs-1

                          i1=n-isent
                          i2=numprocs-islave
                       
                          nloop=i1/i2
                          if(mod(i1,i2).ne.0) nloop=nloop+1

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

                          jsend(01)=ja
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js
                          jsend(17)=jt
                          jsend(18)=ju
                          jsend(19)=jv

                          call mpi_bsend(jsend,19,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

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

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jh),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ji),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jt),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ju),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jv),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          ja=ja+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop
                          jt=jt+nloop
                          ju=ju+nloop
                          jv=jv+nloop

 1017                     isent=isent+nloop   
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,19,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ja=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)
      jt=jsend(17)
      ju=jsend(18)
      jv=jsend(19)

      call mpi_recv(iwkvec(ja),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jh),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ji),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jt),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ju),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jv),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop*3
      kf=ke+nloop*3
      kg=kf+nloop*6
      kh=kg+nloop*10
      ki=kh+nloop
      kj=ki+nloop*3
      kk=kj+ncpods*3
      kl=kk+ncpods*10
      km=kl+ncpods*18

      if(km.gt.nreal) stop 'gradmm ds real overflow'

      call gradovds(nfockdim,ncpods,nloop,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),density,overlap,gradient,wkvec(je),
     &              wkvec(jh),wkvec(jm),wkvec(jn),wkvec(jp),wkvec(jo),
     &              wkvec(jt),wkvec(ju),wkvec(jv),wkvec(ka),wkvec(kb),
     &              wkvec(kc),wkvec(kd),wkvec(ke),wkvec(kf),wkvec(kg),
     &              wkvec(kh),wkvec(ki),wkvec(kj),wkvec(kk),wkvec(kl))

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop*3
      kl=kk+nloop*3
      km=kl+nloop*3
      kn=km+nloop*6
      ko=kn+nloop*6
      kp=ko+nloop*10
      kq=kp+nloop
      kr=kq+nloop*3
      ks=kr+ncpods*3
      kt=ks+ncpods*6
      ku=kt+ncpods*10
      kv=ku+ncpods*18
      kw=kv+ncpods*18

      if(kw.gt.nreal) stop 'gradmm ds real overflow'

      call gradnads(ncpods,nloop,ncenters,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),g,coord,charge,density,gradient,
     &              wkvec(je),wkvec(jh),wkvec(ji),wkvec(jq),wkvec(jr),
     &              wkvec(js),wkvec(jp),wkvec(jt),wkvec(ju),wkvec(jv),
     &              wkvec(jw),wkvec(jx),wkvec(jy),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))

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop*3
      kv=ku+nloop*3
      kw=kv+nloop*3
      kx=kw+nloop*6
      ky=kx+nloop*6
      kz=ky+nloop*10
      k0=kz+nloop
      k1=k0+nloop*3
      k2=k1+ncpods*3
      k3=k2+ncpods*6
      k4=k3+ncpods*10
      k5=k4+ncpods*18
      k6=k5+ncpods*18

      if(k6.gt.nreal) stop 'gradmm ds real overflow'

      call gradcdds(ncpods,nloop,ncds,icdlfunc,icdcfunc,iwkvec(ja),
     &              iwkvec(jb),iwkvec(jc),iwkvec(ia),g,coord,density,
     &              gradient,cdfitc,alphacd,coefscd,wkvec(je),wkvec(jh),
     &              wkvec(ji),wkvec(jq),wkvec(jr),wkvec(js),wkvec(jp),
     &              wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jw),wkvec(jx),
     &              wkvec(jy),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(k0),wkvec(k1),wkvec(k2),
     &              wkvec(k3),wkvec(k4),wkvec(k5))

      ka=jz+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop
      kv=ku+nloop
      kw=kv+nloop
      kx=kw+nloop
      ky=kx+nloop
      kz=ky+nloop
      la=kz+nloop*3
      lb=la+nloop*3
      lc=lb+nloop*3
      ld=lc+nloop*3
      le=ld+nloop*3
      lf=le+nloop*6
      lg=lf+nloop*6
      lh=lg+nloop*6
      li=lh+nloop*6
      lj=li+nloop*10
      lk=lj+nloop*10
      ll=lk+nloop*10
      lm=ll+nloop*3
      ln=lm+nloop*9
      lo=ln+nloop*9
      lp=lo+nloop*18
      lq=lp+nloop*18
      lr=lq+nloop*18
      ls=lr+nloop*36
      lt=ls+nloop*30
      lu=lt+nloop*30
      lv=lu+nloop*60
      lw=lv+nloop
      lx=lw+nloop*10
      ly=lx+nloop*18
      lz=ly+ncpods*3
      ma=lz+ncpods*9
      mb=ma+ncpods*18
      mc=mb+ncpods*6
      md=mc+ncpods*18
      me=md+ncpods*36
      mf=me+ncpods*10
      mg=mf+ncpods*30
      mh=mg+ncpods*60
      mi=mh+ncpods*18
      mj=mi+ncpods*54
      mk=mj+ncpods*108
      ml=mk+ncpods*18

      if(ml.gt.nreal) stop 'gradmm ds real overflow'

      call gradxxds(ncpods,nloop,ncds,ncdspd,icdlfunc,icdcfunc,
     &              iwkvec(ja),iwkvec(jb),iwkvec(jc),iwkvec(ia),g,coord,
     &              density,gradient,cdfitc,alphacd,coefscd,coefpcd,
     &              coefdcd,wkvec(je),wkvec(jh),wkvec(ji),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jp),wkvec(jt),wkvec(ju),
     &              wkvec(jv),wkvec(jw),wkvec(jx),wkvec(jy),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),wkvec(ln),wkvec(lo),
     &              wkvec(lp),wkvec(lq),wkvec(lr),wkvec(ls),wkvec(lt),
     &              wkvec(lu),wkvec(lv),wkvec(lw),wkvec(lx),wkvec(ly),
     &              wkvec(lz),wkvec(ma),wkvec(mb),wkvec(mc),wkvec(md),
     &              wkvec(me),wkvec(mf),wkvec(mg),wkvec(mh),wkvec(mi),
     &              wkvec(mj),wkvec(mk))

      ja=1
      jb=ja+nppodp
      jc=jb+ncpodp*18

      ia=jc+ncpodp*2
      ib=ia+nppodp

      if(ib.gt.ninteger) stop 'gradmm dp integer overflow'

      je=1
      jf=je+nppodp
      jg=jf+nppodp
      jh=jg+nppodp
      ji=jh+nppodp
      jj=ji+nppodp
      jk=jj+nppodp
      jl=jk+nppodp
      jm=jl+nppodp
      jn=jm+nppodp
      jo=jn+nppodp
      jp=jo+nppodp
      jq=jp+nppodp
      jr=jq+nppodp
      js=jr+nppodp
      jt=js+nppodp
      ju=jt+nppodp
      jv=ju+nppodp
      jw=jv+nppodp
      jx=jw+nppodp
      jy=jx+nppodp
      jz=jy+nppodp
      j0=jz+ncpodp
      j1=j0+ncpodp
      j2=j1+ncpodp
      j3=j2+nppodp*16

      if(j2.gt.nreal) stop 'gradmm dp real overflow'

      call setupgdp(ncpodp,n,ncontd,ncontp,nsheld,nshelp,mtloca,iprimvc,
     &              nprimvc,icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),
     &              coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &              wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &              wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &              wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jz),wkvec(j0),
     &              wkvec(j1),wkvec(j2))

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                          call stringcopy(n,wkvec(15*n+1),wkvec(jt))
                          call stringcopy(n,wkvec(16*n+1),wkvec(ju))
                          call stringcopy(n,wkvec(17*n+1),wkvec(jv))
                          call stringcopy(n,wkvec(18*n+1),wkvec(jw))
                          call stringcopy(n,wkvec(19*n+1),wkvec(jx))
                          call stringcopy(n,wkvec(20*n+1),wkvec(jy))
                    endif

      call stringcopy(ncpodp,wkvec(21*n+0*ncpodp+1),wkvec(jz))
      call stringcopy(ncpodp,wkvec(21*n+1*ncpodp+1),wkvec(j0))
      call stringcopy(ncpodp,wkvec(21*n+2*ncpodp+1),wkvec(j1))
 
      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+n
      jy=jx+n
      jz=jy+n
      j0=jz+ncpodp
      j1=j0+ncpodp
      j2=j1+ncpodp

      if(myid.eq.0) then
                          isent=0

                          do 1018 islave=0,numprocs-1

                          i1=n-isent
                          i2=numprocs-islave

                          nloop=i1/i2
                          if(mod(i1,i2).ne.0) nloop=nloop+1

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

                          jsend(01)=ja
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js
                          jsend(17)=jt
                          jsend(18)=ju
                          jsend(19)=jv
                          jsend(20)=jw
                          jsend(21)=jx
                          jsend(22)=jy


                          call mpi_bsend(jsend,22,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

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

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jh),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ji),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jt),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ju),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jv),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jw),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jx),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jy),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          ja=ja+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop
                          jt=jt+nloop
                          ju=ju+nloop
                          jv=jv+nloop
                          jw=jw+nloop
                          jx=jx+nloop
                          jy=jy+nloop

1018                      isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(jsend,22,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ja=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)
      jt=jsend(17)
      ju=jsend(18)
      jv=jsend(19)
      jw=jsend(20)
      jx=jsend(21)
      jy=jsend(22)

      call mpi_recv(iwkvec(ja),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jh),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ji),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jt),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ju),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jv),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jw),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jx),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jy),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop*3
      kf=ke+nloop*3
      kg=kf+nloop*6
      kh=kg+nloop*9
      ki=kh+nloop*10
      kj=ki+nloop*30
      kk=kj+nloop
      kl=kk+nloop*3
      km=kl+nloop*6
      kn=km+ncpodp*9
      ko=kn+ncpodp*30
      kp=ko+ncpodp*54

      if(kp.gt.nreal) stop 'gradmm dp real overflow'

      call gradovdp(nfockdim,ncpodp,nloop,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),density,overlap,gradient,wkvec(je),
     &              wkvec(jh),wkvec(jm),wkvec(jn),wkvec(jp),wkvec(jo),
     &              wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jw),wkvec(jx),
     &              wkvec(jy),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))

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop*3
      km=kl+nloop*3
      kn=km+nloop*3
      ko=kn+nloop*3
      kp=ko+nloop*6
      kq=kp+nloop*6
      kr=kq+nloop*6
      ks=kr+nloop*10
      kt=ks+nloop*10
      ku=kt+nloop*15
      kv=ku+nloop
      kw=kv+nloop*6
      kx=kw+ncpodp*3
      ky=kx+ncpodp*6
      kz=ky+ncpodp*9
      k0=kz+ncpodp*10
      k1=k0+ncpodp*18
      k2=k1+ncpodp*15
      k3=k2+ncpodp*30
      k4=k3+ncpodp*36
      k5=k4+ncpodp*54

      if(k5.gt.nreal) stop 'gradmm dp real overflow'

      call gradnadp(ncpodp,nloop,ncenters,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),g,coord,charge,density,gradient,
     &              wkvec(je),wkvec(jh),wkvec(ji),wkvec(jq),wkvec(jr),
     &              wkvec(js),wkvec(jp),wkvec(jt),wkvec(ju),wkvec(jv),
     &              wkvec(jz),wkvec(j0),wkvec(j1),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(k0),
     &              wkvec(k1),wkvec(k2),wkvec(k3),wkvec(k4))

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop
      kv=ku+nloop*3
      kw=kv+nloop*3
      kx=kw+nloop*3
      ky=kx+nloop*3
      kz=ky+nloop*6
      la=kz+nloop*6
      lb=la+nloop*6
      lc=lb+nloop*10
      ld=lc+nloop*10
      le=ld+nloop*15
      lf=le+nloop
      lg=lf+nloop*6
      lh=lg+ncpodp*3
      li=lh+ncpodp*6
      lj=li+ncpodp*9
      lk=lj+ncpodp*10
      ll=lk+ncpodp*18
      lm=ll+ncpodp*15
      ln=lm+ncpodp*30
      lo=ln+ncpodp*36
      lp=lo+ncpodp*54

      if(lp.gt.nreal) stop 'gradmm dp real overflow'

      call gradcddp(ncpodp,nloop,ncds,icdlfunc,icdcfunc,iwkvec(ja),
     &              iwkvec(jb),iwkvec(jc),iwkvec(ia),g,coord,density,
     &              gradient,cdfitc,alphacd,coefscd,wkvec(je),wkvec(jh),
     &              wkvec(ji),wkvec(jq),wkvec(jr),wkvec(js),wkvec(jp),
     &              wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jz),wkvec(j0),
     &              wkvec(j1),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),
     &              wkvec(ln),wkvec(lo))

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop
      kv=ku+nloop
      kw=kv+nloop
      kx=kw+nloop
      ky=kx+nloop
      kz=ky+nloop
      la=kz+nloop
      lb=la+nloop*3
      lc=lb+nloop*3
      ld=lc+nloop*3
      le=ld+nloop*3
      lf=le+nloop*3
      lg=lf+nloop*3
      lh=lg+nloop*6
      li=lh+nloop*6
      lj=li+nloop*6
      lk=lj+nloop*6
      ll=lk+nloop*6
      lm=ll+nloop*10
      ln=lm+nloop*10
      lo=ln+nloop*10
      lp=lo+nloop*10
      lq=lp+nloop*15
      lr=lq+nloop*15
      ls=lr+nloop*15
      lt=ls+nloop*3
      lu=lt+nloop*9
      lv=lu+nloop*9
      lw=lv+nloop*18
      lx=lw+nloop*18
      ly=lx+nloop*18
      lz=ly+nloop*36
      ma=lz+nloop*30
      mb=ma+nloop*30
      mc=mb+nloop*60
      md=mc+nloop*45
      me=md+nloop*45
      mf=me+nloop*90
      mg=mf+nloop
      mh=mg+nloop*15
      mi=mh+nloop*30
      mj=mi+ncpodp*3
      mk=mj+ncpodp*9
      ml=mk+ncpodp*18
      mm=ml+ncpodp*6
      mn=mm+ncpodp*18
      mo=mn+ncpodp*36
      mp=mo+ncpodp*9
      mq=mp+ncpodp*27
      mr=mq+ncpodp*54
      ms=mr+ncpodp*10
      mt=ms+ncpodp*30
      mu=mt+ncpodp*60
      mv=mu+ncpodp*18
      mw=mv+ncpodp*54
      mx=mw+ncpodp*108
      my=mx+ncpodp*15
      mz=my+ncpodp*45
      na=mz+ncpodp*90
      nb=na+ncpodp*30
      nc=nb+ncpodp*90
      nd=nc+ncpodp*180
      ne=nd+ncpodp*36
      nf=ne+ncpodp*108
      ng=nf+ncpodp*216
      nh=ng+ncpodp*54

      if(nh.gt.nreal) stop 'gradmm dp real overflow'

      call gradxxdp(ncpodp,nloop,ncds,ncdspd,icdlfunc,icdcfunc,
     &              iwkvec(ja),iwkvec(jb),iwkvec(jc),iwkvec(ia),g,coord,
     &              density,gradient,cdfitc,alphacd,coefscd,coefpcd,
     &              coefdcd,wkvec(je),wkvec(jh),wkvec(ji),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jp),wkvec(jt),wkvec(ju),
     &              wkvec(jv),wkvec(jz),wkvec(j0),wkvec(j1),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),wkvec(ln),wkvec(lo),
     &              wkvec(lp),wkvec(lq),wkvec(lr),wkvec(ls),wkvec(lt),
     &              wkvec(lu),wkvec(lv),wkvec(lw),wkvec(lx),wkvec(ly),
     &              wkvec(lz),wkvec(ma),wkvec(mb),wkvec(mc),wkvec(md),
     &              wkvec(me),wkvec(mf),wkvec(mg),wkvec(mh),wkvec(mi),
     &              wkvec(mj),wkvec(mk),wkvec(ml),wkvec(mm),wkvec(mn),
     &              wkvec(mo),wkvec(mp),wkvec(mq),wkvec(mr),wkvec(ms),
     &              wkvec(mt),wkvec(mu),wkvec(mv),wkvec(mw),wkvec(mx),
     &              wkvec(my),wkvec(mz),wkvec(na),wkvec(nb),wkvec(nc),
     &              wkvec(nd),wkvec(ne),wkvec(nf),wkvec(ng))

      ja=1
      jb=ja+ncpodd
      jc=jb+nppodd
      jd=jc+ncpodd*36

      ia=jd+ncpodd*2
      ib=ia+nppodd

      if(ib.gt.ninteger) stop 'gradmm dd integer overflow'

      je=1
      jf=je+nppodd
      jg=jf+nppodd
      jh=jg+nppodd
      ji=jh+nppodd
      jj=ji+nppodd
      jk=jj+nppodd
      jl=jk+nppodd
      jm=jl+nppodd
      jn=jm+nppodd
      jo=jn+nppodd
      jp=jo+nppodd
      jq=jp+nppodd
      jr=jq+nppodd
      js=jr+nppodd
      jt=js+nppodd
      ju=jt+nppodd
      jv=ju+nppodd
      jw=jv+nppodd
      jx=jw+nppodd
      jy=jx+nppodd
      jz=jy+nppodd
      j0=jz+ncpodd
      j1=j0+ncpodd
      j2=j1+ncpodd
      j3=j2+nppodd*16

      if(j3.gt.nreal) stop 'gradmm dd real overflow'

      call setupgdd(ncpodd,n,ncontd,nsheld,mtloca,iprimvc,nprimvc,
     &              icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),iwkvec(jd),
     &              coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &              wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &              wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &              wkvec(jw),wkvec(jx),wkvec(jy),wkvec(jz),wkvec(j0),
     &              wkvec(j1),wkvec(j2))

      if(myid.eq.0) then
                          call stringcopy(n,wkvec(00*n+1),wkvec(je))
                          call stringcopy(n,wkvec(01*n+1),wkvec(jf))
                          call stringcopy(n,wkvec(02*n+1),wkvec(jg))
                          call stringcopy(n,wkvec(03*n+1),wkvec(jh))
                          call stringcopy(n,wkvec(04*n+1),wkvec(ji))
                          call stringcopy(n,wkvec(05*n+1),wkvec(jj))
                          call stringcopy(n,wkvec(06*n+1),wkvec(jk))
                          call stringcopy(n,wkvec(07*n+1),wkvec(jl))
                          call stringcopy(n,wkvec(08*n+1),wkvec(jm))
                          call stringcopy(n,wkvec(09*n+1),wkvec(jn))
                          call stringcopy(n,wkvec(10*n+1),wkvec(jo))
                          call stringcopy(n,wkvec(11*n+1),wkvec(jp))
                          call stringcopy(n,wkvec(12*n+1),wkvec(jq))
                          call stringcopy(n,wkvec(13*n+1),wkvec(jr))
                          call stringcopy(n,wkvec(14*n+1),wkvec(js))
                          call stringcopy(n,wkvec(15*n+1),wkvec(jt))
                          call stringcopy(n,wkvec(16*n+1),wkvec(ju))
                          call stringcopy(n,wkvec(17*n+1),wkvec(jv))
                          call stringcopy(n,wkvec(18*n+1),wkvec(jw))
                          call stringcopy(n,wkvec(19*n+1),wkvec(jx))
                          call stringcopy(n,wkvec(20*n+1),wkvec(jy))
                    endif

      call stringcopy(ncpodd,wkvec(21*n+0*ncpodd+1),wkvec(jz))
      call stringcopy(ncpodd,wkvec(21*n+1*ncpodd+1),wkvec(j0))
      call stringcopy(ncpodd,wkvec(21*n+2*ncpodd+1),wkvec(j1))
 
      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+n
      jy=jx+n
      jz=jy+n
      j0=jz+ncpodd
      j1=j0+ncpodd
      j2=j1+ncpodd

      if(myid.eq.0) then
                          isent=0

                          do 1019 islave=0,numprocs-1

                          i1=n-isent
                          i2=numprocs-islave

                          nloop=i1/i2
                          if(mod(i1,i2).ne.0) nloop=nloop+1

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

                          jsend(01)=jb
                          jsend(02)=je
                          jsend(03)=jf
                          jsend(04)=jg
                          jsend(05)=jh
                          jsend(06)=ji
                          jsend(07)=jj
                          jsend(08)=jk
                          jsend(09)=jl
                          jsend(10)=jm
                          jsend(11)=jn
                          jsend(12)=jo
                          jsend(13)=jp
                          jsend(14)=jq
                          jsend(15)=jr
                          jsend(16)=js
                          jsend(17)=jt
                          jsend(18)=ju
                          jsend(19)=jv
                          jsend(20)=jw
                          jsend(21)=jx
                          jsend(22)=jy


                          call mpi_bsend(jsend,22,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          call mpi_bsend(iwkvec(jb),nloop,
     &                                   mpi_integer,islave,islave,
     &                                   mpi_comm_world,ierr)

                          call mpi_bsend(wkvec(je),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jh),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ji),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jm),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jn),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jo),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jp),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jq),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jr),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(js),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jt),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(ju),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jv),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jw),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jx),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)
                          call mpi_bsend(wkvec(jy),nloop,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

                          jb=jb+nloop

                          je=je+nloop
                          jf=jf+nloop
                          jg=jg+nloop
                          jh=jh+nloop
                          ji=ji+nloop
                          jj=jj+nloop
                          jk=jk+nloop
                          jl=jl+nloop
                          jm=jm+nloop
                          jn=jn+nloop
                          jo=jo+nloop
                          jp=jp+nloop
                          jq=jq+nloop
                          jr=jr+nloop
                          js=js+nloop
                          jt=jt+nloop
                          ju=ju+nloop
                          jv=jv+nloop
                          jw=jw+nloop
                          jx=jx+nloop
                          jy=jy+nloop

1019                      isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierr)

      call mpi_recv(jsend,22,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierr)

      jb=jsend(01)
      je=jsend(02)
      jf=jsend(03)
      jg=jsend(04)
      jh=jsend(05)
      ji=jsend(06)
      jj=jsend(07)
      jk=jsend(08)
      jl=jsend(09)
      jm=jsend(10)
      jn=jsend(11)
      jo=jsend(12)
      jp=jsend(13)
      jq=jsend(14)
      jr=jsend(15)
      js=jsend(16)
      jt=jsend(17)
      ju=jsend(18)
      jv=jsend(19)
      jw=jsend(20)
      jx=jsend(21)
      jy=jsend(22)

      call mpi_recv(iwkvec(jb),nloop,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      call mpi_recv(wkvec(je),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jh),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ji),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jm),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jn),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jo),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jp),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jq),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jr),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(js),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jt),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(ju),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jv),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jw),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jx),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)
      call mpi_recv(wkvec(jy),nloop,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierror)

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop*3
      kf=ke+nloop*3
      kg=kf+nloop*3
      kh=kg+nloop*6
      ki=kh+nloop*6
      kj=ki+nloop*6
      kk=kj+nloop*10
      kl=kk+nloop*10
      km=kl+nloop*18
      kn=km+nloop*18
      ko=kn+nloop*30
      kp=ko+nloop*60
      kq=kp+nloop
      kr=kq+nloop*10
      ks=kr+nloop*18
      kt=ks+ncpodd*18
      ku=kt+ncpodd*60
      kv=ku+ncpodd*108

      if(kv.gt.nreal) stop 'gradmm dd real overflow'

      call gradovdd(nfockdim,ncpodd,nloop,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),iwkvec(jd),density,overlap,gradient,
     &              wkvec(je),wkvec(jh),wkvec(ji),wkvec(jm),wkvec(jn),
     &              wkvec(jp),wkvec(jo),wkvec(jt),wkvec(ju),wkvec(jv),
     &              wkvec(jw),wkvec(jx),wkvec(jy),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))

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop*3
      kn=km+nloop*3
      ko=kn+nloop*3
      kp=ko+nloop*3
      kq=kp+nloop*3
      kr=kq+nloop*6
      ks=kr+nloop*6
      kt=ks+nloop*6
      ku=kt+nloop*6
      kv=ku+nloop*10
      kw=kv+nloop*10
      kx=kw+nloop*10
      ky=kx+nloop*15
      kz=ky+nloop*15
      la=kz+nloop*21
      lb=la+nloop
      lc=lb+nloop*10
      ld=lc+ncpodd*3
      le=ld+ncpodd*6
      lf=le+ncpodd*9
      lg=lf+ncpodd*10
      lh=lg+ncpodd*18
      li=lh+ncpodd*18
      lj=li+ncpodd*15
      lk=lj+ncpodd*30
      ll=lk+ncpodd*36
      lm=ll+ncpodd*21
      ln=lm+ncpodd*45
      lo=ln+ncpodd*60
      lp=lo+ncpodd*60
      lq=lp+ncpodd*108

      if(lq.gt.nreal) stop 'gradmm dd real overflow'

      call gradnadd(ncpodd,nloop,ncenters,iwkvec(ja),iwkvec(jb),
     &              iwkvec(jc),iwkvec(jd),g,coord,charge,density,
     &              gradient,wkvec(je),wkvec(jh),wkvec(ji),wkvec(jq),
     &              wkvec(jr),wkvec(js),wkvec(jp),wkvec(jt),wkvec(ju),
     &              wkvec(jv),wkvec(jz),wkvec(j0),wkvec(j1),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),wkvec(ln),wkvec(lo),
     &              wkvec(lp))

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop
      kv=ku+nloop
      kw=kv+nloop*3
      kx=kw+nloop*3
      ky=kx+nloop*3
      kz=ky+nloop*3
      la=kz+nloop*3
      lb=la+nloop*6
      lc=lb+nloop*6
      ld=lc+nloop*6
      le=ld+nloop*6
      lf=le+nloop*10
      lg=lf+nloop*10
      lh=lg+nloop*10
      li=lh+nloop*15
      lj=li+nloop*15
      lk=lj+nloop*21
      ll=lk+nloop
      lm=ll+nloop*10
      ln=lm+ncpodd*3
      lo=ln+ncpodd*6
      lp=lo+ncpodd*9
      lq=lp+ncpodd*10
      lr=lq+ncpodd*18
      ls=lr+ncpodd*18
      lt=ls+ncpodd*15
      lu=lt+ncpodd*30
      lv=lu+ncpodd*36
      lw=lv+ncpodd*21
      lx=lw+ncpodd*45
      ly=lx+ncpodd*60
      lz=ly+ncpodd*60
      l0=lz+ncpodd*108

      if(l0.gt.nreal) stop 'gradmm dd real overflow'

      call gradcddd(ncpodd,nloop,ncds,icdlfunc,icdcfunc,iwkvec(ja),
     &              iwkvec(jb),iwkvec(jc),iwkvec(jd),iwkvec(ia),g,coord,
     &              density,gradient,cdfitc,alphacd,coefscd,wkvec(je),
     &              wkvec(jh),wkvec(ji),wkvec(jq),wkvec(jr),wkvec(js),
     &              wkvec(jp),wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jz),
     &              wkvec(j0),wkvec(j1),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),wkvec(ln),wkvec(lo),wkvec(lp),wkvec(lq),
     &              wkvec(lr),wkvec(ls),wkvec(lt),wkvec(lu),wkvec(lv),
     &              wkvec(lw),wkvec(lx),wkvec(ly),wkvec(lz))

      ka=j2+1
      kb=ka+nloop
      kc=kb+nloop
      kd=kc+nloop
      ke=kd+nloop
      kf=ke+nloop
      kg=kf+nloop
      kh=kg+nloop
      ki=kh+nloop
      kj=ki+nloop
      kk=kj+nloop
      kl=kk+nloop
      km=kl+nloop
      kn=km+nloop
      ko=kn+nloop
      kp=ko+nloop
      kq=kp+nloop
      kr=kq+nloop
      ks=kr+nloop
      kt=ks+nloop
      ku=kt+nloop
      kv=ku+nloop
      kw=kv+nloop
      kx=kw+nloop
      ky=kx+nloop
      kz=ky+nloop
      la=kz+nloop
      lb=la+nloop
      lc=lb+nloop*3
      ld=lc+nloop*3
      le=ld+nloop*3
      lf=le+nloop*3
      lg=lf+nloop*3
      lh=lg+nloop*3
      li=lh+nloop*3
      lj=li+nloop*6
      lk=lj+nloop*6
      ll=lk+nloop*6
      lm=ll+nloop*6
      ln=lm+nloop*6
      lo=ln+nloop*6
      lp=lo+nloop*10
      lq=lp+nloop*10
      lr=lq+nloop*10
      ls=lr+nloop*10
      lt=ls+nloop*10
      lu=lt+nloop*15
      lv=lu+nloop*15
      lw=lv+nloop*15
      lx=lw+nloop*15
      ly=lx+nloop*21
      lz=ly+nloop*21
      ma=lz+nloop*21
      mb=ma+nloop*3
      mc=mb+nloop*9
      md=mc+nloop*9
      me=md+nloop*18
      mf=me+nloop*18
      mg=mf+nloop*18
      mh=mg+nloop*36
      mi=mh+nloop*30
      mj=mi+nloop*30
      mk=mj+nloop*60
      ml=mk+nloop*45
      mm=ml+nloop*45
      mn=mm+nloop*90
      mo=mn+nloop*63
      mp=mo+nloop*63
      mq=mp+nloop*126
      mr=mq+nloop
      ms=mr+nloop*21
      mt=ms+nloop*45
      mu=mt+ncpodd*3
      mv=mu+ncpodd*9
      mw=mv+ncpodd*18
      mx=mw+ncpodd*6
      my=mx+ncpodd*18
      mz=my+ncpodd*36
      na=mz+ncpodd*9
      nb=na+ncpodd*27
      nc=nb+ncpodd*54
      nd=nc+ncpodd*10
      ne=nd+ncpodd*30
      nf=ne+ncpodd*60
      ng=nf+ncpodd*18
      nh=ng+ncpodd*54
      ni=nh+ncpodd*108
      nj=ni+ncpodd*18
      nk=nj+ncpodd*54
      nl=nk+ncpodd*108
      nm=nl+ncpodd*15
      nn=nm+ncpodd*45
      no=nn+ncpodd*90
      np=no+ncpodd*30
      nq=np+ncpodd*90
      nr=nq+ncpodd*180
      ns=nr+ncpodd*36
      nt=ns+ncpodd*108
      nu=nt+ncpodd*216
      nv=nu+ncpodd*21
      nw=nv+ncpodd*63
      nx=nw+ncpodd*126
      ny=nx+ncpodd*45
      nz=ny+ncpodd*135
      n0=nz+ncpodd*270
      n1=n0+ncpodd*60
      n2=n1+ncpodd*180
      n3=n2+ncpodd*360
      n4=n3+ncpodd*60
      n5=n4+ncpodd*180
      n6=n5+ncpodd*360
      n7=n6+ncpodd*108

      if(n7.gt.nreal) stop 'gradmm dd real overflow'

      call gradxxdd(ncpodd,nloop,ncds,ncdspd,icdlfunc,icdcfunc,
     &              iwkvec(ja),iwkvec(jb),iwkvec(jc),iwkvec(jd),
     &              iwkvec(ia),g,coord,density,gradient,cdfitc,alphacd,
     &              coefscd,coefpcd,coefdcd,wkvec(je),wkvec(jh),
     &              wkvec(ji),wkvec(jq),wkvec(jr),wkvec(js),wkvec(jp),
     &              wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jz),wkvec(j0),
     &              wkvec(j1),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),
     &              wkvec(ln),wkvec(lo),wkvec(lp),wkvec(lq),wkvec(lr),
     &              wkvec(ls),wkvec(lt),wkvec(lu),wkvec(lv),wkvec(lw),
     &              wkvec(lx),wkvec(ly),wkvec(lz),wkvec(ma),wkvec(mb),
     &              wkvec(mc),wkvec(md),wkvec(me),wkvec(mf),wkvec(mg),
     &              wkvec(mh),wkvec(mi),wkvec(mj),wkvec(mk),wkvec(ml),
     &              wkvec(mm),wkvec(mn),wkvec(mo),wkvec(mp),wkvec(mq),
     &              wkvec(mr),wkvec(ms),wkvec(mt),wkvec(mu),wkvec(mv),
     &              wkvec(mw),wkvec(mx),wkvec(my),wkvec(mz),wkvec(na),
     &              wkvec(nb),wkvec(nc),wkvec(nd),wkvec(ne),wkvec(nf),
     &              wkvec(ng),wkvec(nh),wkvec(ni),wkvec(nj),wkvec(nk),
     &              wkvec(nl),wkvec(nm),wkvec(nn),wkvec(no),wkvec(np),
     &              wkvec(nq),wkvec(nr),wkvec(ns),wkvec(nt),wkvec(nu),
     &              wkvec(nv),wkvec(nw),wkvec(nx),wkvec(ny),wkvec(nz),
     &              wkvec(n0),wkvec(n1),wkvec(n2),wkvec(n3),wkvec(n4),
     &              wkvec(n5),wkvec(n6))

 2001 continue

      call normo(nfockdim,nconts,ncontp,ncontd,
     &           nshels,nshelp,nsheld,mtloca,overlap)

      call mpi_reduce(overlap(1,1),wkvec,nfockdim,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierror)

      if(myid.eq.0) call stringcopy(nfockdim,overlap(1,1),wkvec)

      call mpi_reduce(overlap(1,2),wkvec,nfockdim,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierror)

      if(myid.eq.0) call stringcopy(nfockdim,overlap(1,2),wkvec)

      call mpi_reduce(overlap(1,3),wkvec,nfockdim,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierror)

      if(myid.eq.0) call stringcopy(nfockdim,overlap(1,3),wkvec)

      call mpi_reduce(gradient,wkvec,3*ncenters,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)

      if(myid.eq.0) call stringcopy(3*ncenters,gradient,wkvec)

      call mpi_bcast(overlap,3*nfockdim,mpi_real8,0,mpi_comm_world,ierr)

      return
      end
