      subroutine cdgrxs(ncds,ncdspd,ilfunc,icfunc,cdfitc,alpha,coeffs,
     &                  coeffp,coeffd,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(4)
 
      dimension ilfunc(*),icfunc(*)

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

      data one,two,three/1.0d0,2.0d0,3.0d0/
      data four,five,forty9/4.0d0,5.0d0,49.0d0/
      data zero,pi/0.0d0,3.1415926535898d0/

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

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

      ncombos=ncds*ncdspd

      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,ncdspd
      do 1001 j=1,ncds

      ij=ij+1

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

      ilabel=ilfunc(ncds+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(4,u,g)

      s0=g(1)*const*sqrt(zetaab)
      s1=g(2)*const*sqrt(zetaab)
      s2=g(3)*const*sqrt(zetaab)
      s3=g(4)*const*sqrt(zetaab)

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

      px1=wxpx*s2
      py1=wypy*s2
      pz1=wzpz*s2

      px2=wxpx*s3
      py2=wypy*s3
      pz2=wzpz*s3

      cfact=(s0-zetab*zetaab*s1)/(two*zetaa)

      dxx0=wxpx*px1+cfact
      dxy0=wypy*px1
      dxz0=wzpz*px1
      dyy0=wypy*py1+cfact
      dyz0=wzpz*py1
      dzz0=wzpz*pz1+cfact

      cfact=(s1-zetab*zetaab*s2)/(two*zetaa)

      dxx1=wxpx*px2+cfact
      dxy1=wypy*px2
      dxz1=wzpz*px2
      dyy1=wypy*py2+cfact
      dyz1=wzpz*py2
      dzz1=wzpz*pz2+cfact

      cfactx=(px0-zetab*zetaab*px1)/(two*zetaa)
      cfacty=(py0-zetab*zetaab*py1)/(two*zetaa)
      cfactz=(pz0-zetab*zetaab*pz1)/(two*zetaa)

      fxxx0=wxpx*dxx1+cfactx+cfactx
      fxxy0=wypy*dxx1
      fxxz0=wzpz*dxx1
      fxyy0=wypy*dxy1+cfactx
      fxyz0=wzpz*dxy1
      fxzz0=wzpz*dxz1+cfactx
      fyyy0=wypy*dyy1+cfacty+cfacty
      fyyz0=wzpz*dyy1
      fyzz0=wzpz*dyz1+cfacty
      fzzz0=wzpz*dzz1+cfactz+cfactz

      n=ncds+(i-1)*10

      cfa01=cdfitc(n+01)
      cfa02=cdfitc(n+02)
      cfa03=cdfitc(n+03)
      cfa04=cdfitc(n+04)
      cfa05=cdfitc(n+05)/sqrt493
      cfa06=cdfitc(n+06)
      cfa07=cdfitc(n+07)
      cfa08=cdfitc(n+08)/sqrt493
      cfa09=cdfitc(n+09)
      cfa10=cdfitc(n+10)/sqrt493

      cfb01=cdfitc(j)

      css=cfa01*cfb01*coeffs(ilabel)*coeffs(jlabel)
      cpx=cfa02*cfb01*coeffp(ilabel)*coeffs(jlabel)
      cpy=cfa03*cfb01*coeffp(ilabel)*coeffs(jlabel)
      cpz=cfa04*cfb01*coeffp(ilabel)*coeffs(jlabel)
      cxx=cfa05*cfb01*coeffd(ilabel)*coeffs(jlabel)
      cxy=cfa06*cfb01*coeffd(ilabel)*coeffs(jlabel)
      cxz=cfa07*cfb01*coeffd(ilabel)*coeffs(jlabel)
      cyy=cfa08*cfb01*coeffd(ilabel)*coeffs(jlabel)
      cyz=cfa09*cfb01*coeffd(ilabel)*coeffs(jlabel)
      czz=cfa10*cfb01*coeffd(ilabel)*coeffs(jlabel)

      factorx=zero
      factory=zero
      factorz=zero

      factorx=factorx-two*zetaa*(css*px0
     &                          +cpx*dxx0  +cpy*dxy0  +cpz*dxz0
     &                          +cxx*fxxx0 +cxy*fxxy0 +cxz*fxxz0
     &                          +cyy*fxyy0 +cyz*fxyz0 +czz*fxzz0)

      factory=factory-two*zetaa*(css*py0
     &                          +cpx*dxy0  +cpy*dyy0  +cpz*dyz0
     &                          +cxx*fxxy0 +cxy*fxyy0 +cxz*fxyz0
     &                          +cyy*fyyy0 +cyz*fyyz0 +czz*fyzz0)

      factorz=factorz-two*zetaa*(css*pz0
     &                          +cpx*dxz0  +cpy*dyz0  +cpz*dzz0
     &                          +cxx*fxxz0 +cxy*fxyz0 +cxz*fxzz0
     &                          +cyy*fyyz0 +cyz*fyzz0 +czz*fzzz0)

      factorx=factorx+cpx*s0+cxx*px0+cxy*py0+cxz*pz0+cxx*px0
      factory=factory+cpy*s0+cxy*px0+cyy*py0+cyz*pz0+cyy*py0
      factorz=factorz+cpz*s0+cxz*px0+cyz*py0+czz*pz0+czz*pz0

      gradient(1,nc1)=gradient(1,nc1)+factorx
      gradient(2,nc1)=gradient(2,nc1)+factory
      gradient(3,nc1)=gradient(3,nc1)+factorz

      gradient(1,nc2)=gradient(1,nc2)-factorx
      gradient(2,nc2)=gradient(2,nc2)-factory
      gradient(3,nc2)=gradient(3,nc2)-factorz

 1001 continue

      return
      end
