      subroutine gradcdss(m,n,ncds,ilfunc,icfunc,icontr,inddyy,icentr,
     &                    iicontr,g,coord,density,gradient,cdfitc,
     &                    alphaaux,coeffs,zeta,twozetaa,twozetab,px,py,
     &                    pz,base,pxax,pyay,pzaz,axbx,ayby,azbz,zzeta,
     &                    ttwozetaa,ttwozetab,zetarec,ppx,ppy,ppz,ppxax,
     &                    ppyay,ppzaz,wxpx,wypy,wzpz,u,sss0,sss1,pss0,
     &                    sscond,pscond,spcond,grad)

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)

      dimension ilfunc(*),icfunc(*),icontr(*),inddyy(*),icentr(m,*),
     &          iicontr(*)

      dimension g(*),coord(3,*),density(*),gradient(3,*),cdfitc(*),
     &          alphaaux(*),coeffs(*),zeta(*),twozetaa(*),twozetab(*),
     &          px(*),py(*),pz(*),base(*),pxax(*),pyay(*),pzaz(*),
     &          axbx(*),ayby(*),azbz(*),zzeta(*),ttwozetaa(*),
     &          ttwozetab(*),zetarec(*),ppx(*),ppy(*),ppz(*),ppxax(*),
     &          ppyay(*),ppzaz(*),wxpx(*),wypy(*),wzpz(*),u(*),sss0(*),
     &          sss1(*),pss0(n,*),sscond(*),pscond(m,*),spcond(m,*),
     &          grad(m,*)

      data one,two,three,four,five,seven/1.0d0,2.0d0,3.0d0,
     &                                   4.0d0,5.0d0,7.0d0/
      data zero,half,thrhlf,pi/0.0d0,0.5d0,1.5d0,3.1415926535898d0/
      data calcut/1.0d-16/

      twoeconst=two*pi**(five/two)

      do 1001 naux=1,ncds

      ilaux=ilfunc(naux)

      nc=icfunc(ilaux)

      ccx=coord(1,nc)
      ccy=coord(2,nc)
      ccz=coord(3,nc)

      zetaux=alphaaux(ilaux)

      cdfits=cdfitc(naux)*coeffs(ilaux)

      do 1002 i=1,n
 1002 zetarec(i)=one/(zeta(i)+zetaux)

      do 1003 i=1,n
      pxcx=px(i)-ccx
      pycy=py(i)-ccy
      pzcz=pz(i)-ccz
 1003 u(i)=zeta(i)*zetaux*zetarec(i)*(pxcx*pxcx+pycy*pycy+pzcz*pzcz)

      do 1004 i=1,n
      call augg(2,u(i),g)
      sss0(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(1)/zeta(i)/zetaux
 1004 sss1(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(2)/zeta(i)/zetaux

      nn=0

      do 1005 i=1,n
      if(abs(sss0(i)).gt.calcut) then
                                       nn=nn+1

                                       iicontr(nn)=icontr(i)

                                       zzeta(nn)=zeta(i)

                                       ttwozetaa(nn)=twozetaa(i)
                                       ttwozetab(nn)=twozetab(i)

                                       ppx(nn)=px(i)
                                       ppy(nn)=py(i)
                                       ppz(nn)=pz(i)

                                       ppxax(nn)=pxax(i)
                                       ppyay(nn)=pyay(i)
                                       ppzaz(nn)=pzaz(i)

                                       zetarec(nn)=zetarec(i)

                                       sss0(nn)=sss0(i)
                                       sss1(nn)=sss1(i)
                                 endif
 1005 continue

      do 1006 i=1,nn
      wxpx(i)=((zzeta(i)*ppx(i)+zetaux*ccx)*zetarec(i))-ppx(i)
      wypy(i)=((zzeta(i)*ppy(i)+zetaux*ccy)*zetarec(i))-ppy(i)
 1006 wzpz(i)=((zzeta(i)*ppz(i)+zetaux*ccz)*zetarec(i))-ppz(i)

      do 1007 i=1,nn
      pss0(i,1)=ppxax(i)*sss0(i)+wxpx(i)*sss1(i)
      pss0(i,2)=ppyay(i)*sss0(i)+wypy(i)*sss1(i)
 1007 pss0(i,3)=ppzaz(i)*sss0(i)+wzpz(i)*sss1(i)

      do 1008 j=1,3
      do 1008 i=1,m
 1008 pscond(i,j)=zero

      do 1009 j=1,3
      do 1009 i=1,nn
 1009 pscond(iicontr(i),j)=pscond(iicontr(i),j)+pss0(i,j)*ttwozetaa(i)

      do 1010 j=1,3
      do 1010 i=1,m
 1010 grad(i,j)=pscond(i,j)*density(inddyy(i))*cdfits

      do 1011 j=1,3
      do 1011 i=1,m
 1011 gradient(j,nc)=gradient(j,nc)-grad(i,j)

      do 1012 j=1,3
      do 1012 i=1,m
 1012 gradient(j,icentr(i,1))=gradient(j,icentr(i,1))+grad(i,j)

      do 1013 i=1,m
 1013 sscond(i)=zero

      do 1014 j=1,3
      do 1014 i=1,m
 1014 pscond(i,j)=zero

      do 1015 i=1,nn
 1015 sscond(iicontr(i))=sscond(iicontr(i))+sss0(i)*ttwozetab(i)

      do 1016 j=1,3
      do 1016 i=1,nn
 1016 pscond(iicontr(i),j)=pscond(iicontr(i),j)+pss0(i,j)*ttwozetab(i)

      do 1017 i=1,m
      spcond(i,1)=axbx(i)*sscond(i)+pscond(i,1)
      spcond(i,2)=ayby(i)*sscond(i)+pscond(i,2)
 1017 spcond(i,3)=azbz(i)*sscond(i)+pscond(i,3)

      do 1018 j=1,3
      do 1018 i=1,m
 1018 grad(i,j)=spcond(i,j)*density(inddyy(i))*cdfits

      do 1019 j=1,3
      do 1019 i=1,m
 1019 gradient(j,nc)=gradient(j,nc)-grad(i,j)

      do 1020 j=1,3
      do 1020 i=1,m
 1020 gradient(j,icentr(i,2))=gradient(j,icentr(i,2))+grad(i,j)

 1001 continue

      return
      end
