      subroutine gradxxss(m,n,ncds,ncdspd,ilfunc,icfunc,icontr,inddyy,
     &                    icentr,iicontr,g,coord,density,gradient,
     &                    cdfitc,alphaaux,coeffs,coeffp,coeffd,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,wxcx,wycy,wzcz,zetac1,zetac2,u,sss0,
     &                    sss1,sss2,sss3,pss0,pss1,pss2,ssp0,ssp1,ssd0,
     &                    psp0,psp1,psd0,c,c1,c2,ssscond,sspcond,
     &                    ssdcond,psscond,pspcond,psdcond,spscond,
     &                    sppcond,spdcond,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 cdfitp(3),cdfitd(6)

      dimension g(*),coord(3,*),density(*),gradient(3,*),cdfitc(*),
     &          alphaaux(*),coeffs(*),coeffp(*),coeffd(*),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(*),wxcx(*),wycy(*),wzcz(*),zetac1(*),
     &          zetac2(*),u(*),sss0(*),sss1(*),sss2(*),sss3(*),
     &          pss0(n,*),pss1(n,*),pss2(n,*),ssp0(n,*),ssp1(n,*),
     &          ssd0(n,*),psp0(n,*),psp1(n,*),psd0(n,*),c(*),c1(n,*),
     &          c2(n,*),ssscond(*),sspcond(m,*),ssdcond(m,*),
     &          psscond(m,*),pspcond(m,*),psdcond(m,*),spscond(m,*),
     &          sppcond(m,*),spdcond(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)

      dnormcd=seven/sqrt(three)

      do 1001 naux=1,ncdspd

      ilaux=ilfunc(ncds+naux)

      nc=icfunc(ilaux)

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

      zetaux=alphaaux(ilaux)

      cdfits=cdfitc(ncds+(naux-1)*10+1)*coeffs(ilaux)

      cdfitp(1)=cdfitc(ncds+(naux-1)*10+2)*coeffp(ilaux)
      cdfitp(2)=cdfitc(ncds+(naux-1)*10+3)*coeffp(ilaux)
      cdfitp(3)=cdfitc(ncds+(naux-1)*10+4)*coeffp(ilaux)

      cdfitd(1)=cdfitc(ncds+(naux-1)*10+05)*coeffd(ilaux)/dnormcd
      cdfitd(2)=cdfitc(ncds+(naux-1)*10+06)*coeffd(ilaux)
      cdfitd(3)=cdfitc(ncds+(naux-1)*10+07)*coeffd(ilaux)
      cdfitd(4)=cdfitc(ncds+(naux-1)*10+08)*coeffd(ilaux)/dnormcd
      cdfitd(5)=cdfitc(ncds+(naux-1)*10+09)*coeffd(ilaux)
      cdfitd(6)=cdfitc(ncds+(naux-1)*10+10)*coeffd(ilaux)/dnormcd

      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(4,u(i),g)
      sss0(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(1)/zeta(i)/zetaux
      sss1(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(2)/zeta(i)/zetaux
      sss2(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(3)/zeta(i)/zetaux
 1004 sss3(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(4)/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)
                                       sss2(nn)=sss2(i)
                                       sss3(nn)=sss3(i)
                                 endif
 1005 continue

      do 1006 i=1,nn
      zetac1(i)=zetaux*zetarec(i)
 1006 zetac2(i)=zzeta(i)*zetarec(i)

      do 1007 i=1,nn
      wxcx(i)=((zzeta(i)*ppx(i)+zetaux*ccx)*zetarec(i))-ccx
      wycy(i)=((zzeta(i)*ppy(i)+zetaux*ccy)*zetarec(i))-ccy
 1007 wzcz(i)=((zzeta(i)*ppz(i)+zetaux*ccz)*zetarec(i))-ccz

      do 1008 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)
 1008 wzpz(i)=((zzeta(i)*ppz(i)+zetaux*ccz)*zetarec(i))-ppz(i)

      do 1009 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)
 1009 pss0(i,3)=ppzaz(i)*sss0(i)+wzpz(i)*sss1(i)

      do 1010 i=1,nn
      pss1(i,1)=ppxax(i)*sss1(i)+wxpx(i)*sss2(i)
      pss1(i,2)=ppyay(i)*sss1(i)+wypy(i)*sss2(i)
 1010 pss1(i,3)=ppzaz(i)*sss1(i)+wzpz(i)*sss2(i)

      do 1011 i=1,nn
      pss2(i,1)=ppxax(i)*sss2(i)+wxpx(i)*sss3(i)
      pss2(i,2)=ppyay(i)*sss2(i)+wypy(i)*sss3(i)
 1011 pss2(i,3)=ppzaz(i)*sss2(i)+wzpz(i)*sss3(i)

      do 1012 i=1,nn
      ssp0(i,1)=wxcx(i)*sss1(i)
      ssp0(i,2)=wycy(i)*sss1(i)
 1012 ssp0(i,3)=wzcz(i)*sss1(i)

      do 1013 i=1,nn
      ssp1(i,1)=wxcx(i)*sss2(i)
      ssp1(i,2)=wycy(i)*sss2(i)
 1013 ssp1(i,3)=wzcz(i)*sss2(i)

      do 1014 i=1,nn
 1014 c(i)=(sss0(i)-zetac2(i)*sss1(i))/(two*zetaux)

      do 1015 i=1,nn
      ssd0(i,1)=wxcx(i)*ssp1(i,1)+c(i)
      ssd0(i,2)=wycy(i)*ssp1(i,1)
      ssd0(i,3)=wzcz(i)*ssp1(i,1)
      ssd0(i,4)=wycy(i)*ssp1(i,2)+c(i)
      ssd0(i,5)=wzcz(i)*ssp1(i,2)
 1015 ssd0(i,6)=wzcz(i)*ssp1(i,3)+c(i)

      do 1016 i=1,nn
 1016 c(i)=half*sss1(i)*zetarec(i)

      do 1017 i=1,nn
      psp0(i,1)=wxcx(i)*pss1(i,1)+c(i)
      psp0(i,2)=wycy(i)*pss1(i,1)
      psp0(i,3)=wzcz(i)*pss1(i,1)
      psp0(i,4)=wxcx(i)*pss1(i,2)
      psp0(i,5)=wycy(i)*pss1(i,2)+c(i)
      psp0(i,6)=wzcz(i)*pss1(i,2)
      psp0(i,7)=wxcx(i)*pss1(i,3)
      psp0(i,8)=wycy(i)*pss1(i,3)
 1017 psp0(i,9)=wzcz(i)*pss1(i,3)+c(i)

      do 1018 i=1,nn
 1018 c(i)=half*sss2(i)*zetarec(i)

      do 1019 i=1,nn
      psp1(i,1)=wxcx(i)*pss2(i,1)+c(i)
      psp1(i,2)=wycy(i)*pss2(i,1)
      psp1(i,3)=wzcz(i)*pss2(i,1)
      psp1(i,4)=wxcx(i)*pss2(i,2)
      psp1(i,5)=wycy(i)*pss2(i,2)+c(i)
      psp1(i,6)=wzcz(i)*pss2(i,2)
      psp1(i,7)=wxcx(i)*pss2(i,3)
      psp1(i,8)=wycy(i)*pss2(i,3)
 1019 psp1(i,9)=wzcz(i)*pss2(i,3)+c(i)

      do 1020 j=1,3
      do 1020 i=1,nn
 1020 c1(i,j)=(pss0(i,j)-zetac2(i)*pss1(i,j))/(two*zetaux)

      do 1021 j=1,3
      do 1021 i=1,nn
 1021 c2(i,j)=half*ssp1(i,j)*zetarec(i)

      do 1022 i=1,nn
      psd0(i,01)=wxcx(i)*psp1(i,1)+c1(i,1)+c2(i,1)
      psd0(i,02)=wycy(i)*psp1(i,1)
      psd0(i,03)=wzcz(i)*psp1(i,1)
      psd0(i,04)=wycy(i)*psp1(i,2)+c1(i,1)
      psd0(i,05)=wzcz(i)*psp1(i,2)
      psd0(i,06)=wzcz(i)*psp1(i,3)+c1(i,1)
      psd0(i,07)=wxcx(i)*psp1(i,4)+c1(i,2)
      psd0(i,08)=wycy(i)*psp1(i,4)        +c2(i,1)
      psd0(i,09)=wzcz(i)*psp1(i,4)
      psd0(i,10)=wycy(i)*psp1(i,5)+c1(i,2)+c2(i,2)
      psd0(i,11)=wzcz(i)*psp1(i,5)
      psd0(i,12)=wzcz(i)*psp1(i,6)+c1(i,2)
      psd0(i,13)=wxcx(i)*psp1(i,7)+c1(i,3)
      psd0(i,14)=wycy(i)*psp1(i,7)
      psd0(i,15)=wzcz(i)*psp1(i,7)        +c2(i,1)
      psd0(i,16)=wycy(i)*psp1(i,8)+c1(i,3)
      psd0(i,17)=wzcz(i)*psp1(i,8)        +c2(i,2)
 1022 psd0(i,18)=wzcz(i)*psp1(i,9)+c1(i,3)+c2(i,3)

      do 1023 j=1,3
      do 1023 i=1,m
 1023 psscond(i,j)=zero

      do 1024 j=1,3
      do 1024 i=1,nn
 1024 psscond(iicontr(i),j)=psscond(iicontr(i),j)
     &                     +pss0(i,j)*ttwozetaa(i)*cdfits

      do 1025 k=1,3
      do 1025 j=1,3
      do 1025 i=1,nn
 1025 psscond(iicontr(i),j)=psscond(iicontr(i),j)
     &                     +psp0(i,3*(j-1)+k)*ttwozetaa(i)*cdfitp(k)

      do 1026 k=1,6
      do 1026 j=1,3
      do 1026 i=1,nn
 1026 psscond(iicontr(i),j)=psscond(iicontr(i),j)
     &                     +psd0(i,6*(j-1)+k)*ttwozetaa(i)*cdfitd(k)

      do 1027 i=1,m
      grad(i,1)=psscond(i,1)*density(inddyy(i))
      grad(i,2)=psscond(i,2)*density(inddyy(i))
 1027 grad(i,3)=psscond(i,3)*density(inddyy(i))

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

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

      do 1030 i=1,m
 1030 ssscond(i)=zero

      do 1031 j=1,3
      do 1031 i=1,m
 1031 psscond(i,j)=zero

      do 1032 i=1,nn
 1032 ssscond(iicontr(i))=ssscond(iicontr(i))
     &                   +sss0(i)*ttwozetab(i)*cdfits

      do 1033 j=1,3
      do 1033 i=1,nn
 1033 ssscond(iicontr(i))=ssscond(iicontr(i))
     &                   +ssp0(i,j)*ttwozetab(i)*cdfitp(j)

      do 1034 j=1,6
      do 1034 i=1,nn
 1034 ssscond(iicontr(i))=ssscond(iicontr(i))
     &                   +ssd0(i,j)*ttwozetab(i)*cdfitd(j)

      do 1035 j=1,3
      do 1035 i=1,nn
 1035 psscond(iicontr(i),j)=psscond(iicontr(i),j)
     &                     +pss0(i,j)*ttwozetab(i)*cdfits

      do 1036 k=1,3
      do 1036 j=1,3
      do 1036 i=1,nn
 1036 psscond(iicontr(i),j)=psscond(iicontr(i),j)
     &                     +psp0(i,3*(j-1)+k)*ttwozetab(i)*cdfitp(k)

      do 1037 k=1,6
      do 1037 j=1,3
      do 1037 i=1,nn
 1037 psscond(iicontr(i),j)=psscond(iicontr(i),j)
     &                     +psd0(i,6*(j-1)+k)*ttwozetab(i)*cdfitd(k)

      do 1038 i=1,m
      spscond(i,1)=axbx(i)*ssscond(i)+psscond(i,1)
      spscond(i,2)=ayby(i)*ssscond(i)+psscond(i,2)
 1038 spscond(i,3)=azbz(i)*ssscond(i)+psscond(i,3)

      do 1039 i=1,m
      grad(i,1)=spscond(i,1)*density(inddyy(i))
      grad(i,2)=spscond(i,2)*density(inddyy(i))
 1039 grad(i,3)=spscond(i,3)*density(inddyy(i))

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

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

 1001 continue

      return
      end
