      subroutine cdgrss(ncds,ilfunc,icfunc,cdfitc,alpha,coeffs,coord,
     &                  gradient)

c  this subroutine written by alain st-amant of the department of
c  pharmaceutical chemistry, university of california, san francisco.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      include "mpif.h"

      dimension istat(mpi_status_size)

      dimension g(2)
 
      dimension ilfunc(*),icfunc(*)

      dimension cdfitc(*),alpha(*),coeffs(*),coord(3,*),gradient(3,*) 

      data one,two,pi/1.0d0,2.0d0,3.1415926535898d0/
      data four,five/4.0d0,5.0d0/

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

      factor=(pi**(five/four)*sqrt(two))**2

      ncombos=(ncds*(ncds+1))/2

      nloop=ncombos/numprocs
      if((nloop*numprocs).lt.ncombos) nloop=nloop+1

      istart=myid*nloop+1
      iend=(myid+1)*nloop

      ij=0

      do 1001 i=1,ncds
      do 1001 j=1,i

      ij=ij+1

      if((ij.lt.istart).or.(ij.gt.iend)) goto 1001

      ilabel=ilfunc(i)
      jlabel=ilfunc(j)

      nc1=icfunc(ilabel)
      nc2=icfunc(jlabel)

      zetaa=alpha(ilabel)
      zetab=alpha(jlabel)

      axbx=coord(1,nc1)-coord(1,nc2)
      ayby=coord(2,nc1)-coord(2,nc2)
      azbz=coord(3,nc1)-coord(3,nc2)

      const=factor/(zetaa*zetab)

      zetaab=one/(zetaa+zetab)

      wxpx=(zetaa*coord(1,nc1)+zetab*coord(1,nc2))*zetaab-coord(1,nc1)
      wypy=(zetaa*coord(2,nc1)+zetab*coord(2,nc2))*zetaab-coord(2,nc1)
      wzpz=(zetaa*coord(3,nc1)+zetab*coord(3,nc2))*zetaab-coord(3,nc1)

      u=zetaa*zetab*zetaab*(axbx*axbx+ayby*ayby+azbz*azbz)

      call augg(2,u,g)

      s1=g(2)*const*sqrt(zetaab)

      px0=wxpx*s1
      py0=wypy*s1
      pz0=wzpz*s1

      cfact=zetaa*cdfitc(i)*cdfitc(j)*coeffs(ilabel)*coeffs(jlabel)

      if(i.ne.j) cfact=cfact+cfact

      gradient(1,nc1)=gradient(1,nc1)-cfact*px0
      gradient(2,nc1)=gradient(2,nc1)-cfact*py0
      gradient(3,nc1)=gradient(3,nc1)-cfact*pz0

      gradient(1,nc2)=gradient(1,nc2)+cfact*px0
      gradient(2,nc2)=gradient(2,nc2)+cfact*py0
      gradient(3,nc2)=gradient(3,nc2)+cfact*pz0

 1001 continue

      return
      end
