      subroutine gradxxds(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,sss4,sss5,pss0,pss1,pss2,pss3,
     &                    pss4,dss0,dss1,dss2,dss3,fss0,fss1,fss2,ssp1,
     &                    psp0,psp1,psd0,dsp0,dsp1,dsd0,fsp0,fsp1,fsd0,
     &                    c,c1,c2,psscond,pspcond,psdcond,dsscond,
     &                    dspcond,dsdcond,fsscond,fspcond,fsdcond,
     &                    dpscond,dppcond,dpdcond,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(m,*),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(*),sss4(*),
     &          sss5(*),pss0(n,*),pss1(n,*),pss2(n,*),pss3(n,*),
     &          pss4(n,*),dss0(n,*),dss1(n,*),dss2(n,*),dss3(n,*),
     &          fss0(n,*),fss1(n,*),fss2(n,*),ssp1(n,*),psp0(n,*),
     &          psp1(n,*),psd0(n,*),dsp0(n,*),dsp1(n,*),dsd0(n,*),
     &          fsp0(n,*),fsp1(n,*),fsd0(n,*),c(*),c1(n,*),c2(n,*),
     &          psscond(m,*),pspcond(m,*),psdcond(m,*),dsscond(m,*),
     &          dspcond(m,*),dsdcond(m,*),fsscond(m,*),fspcond(m,*),
     &          fsdcond(m,*),dpscond(m,*),dppcond(m,*),dpdcond(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(6,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
      sss3(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(4)/zeta(i)/zetaux
      sss4(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(5)/zeta(i)/zetaux
 1004 sss5(i)=base(i)*twoeconst*sqrt(zetarec(i))*g(6)/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)
                                       sss4(nn)=sss4(i)
                                       sss5(nn)=sss5(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
      pss3(i,1)=ppxax(i)*sss3(i)+wxpx(i)*sss4(i)
      pss3(i,2)=ppyay(i)*sss3(i)+wypy(i)*sss4(i)
 1012 pss3(i,3)=ppzaz(i)*sss3(i)+wzpz(i)*sss4(i)

      do 1013 i=1,nn
      pss4(i,1)=ppxax(i)*sss4(i)+wxpx(i)*sss5(i)
      pss4(i,2)=ppyay(i)*sss4(i)+wypy(i)*sss5(i)
 1013 pss4(i,3)=ppzaz(i)*sss4(i)+wzpz(i)*sss5(i)

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

      do 1015 i=1,nn
      dss0(i,1)=ppxax(i)*pss0(i,1)+wxpx(i)*pss1(i,1)+c(i)
      dss0(i,2)=ppyay(i)*pss0(i,1)+wypy(i)*pss1(i,1)
      dss0(i,3)=ppzaz(i)*pss0(i,1)+wzpz(i)*pss1(i,1)
      dss0(i,4)=ppyay(i)*pss0(i,2)+wypy(i)*pss1(i,2)+c(i)
      dss0(i,5)=ppzaz(i)*pss0(i,2)+wzpz(i)*pss1(i,2)
 1015 dss0(i,6)=ppzaz(i)*pss0(i,3)+wzpz(i)*pss1(i,3)+c(i)

      do 1016 i=1,nn
 1016 c(i)=(sss1(i)-zetac1(i)*sss2(i))/(two*zzeta(i))

      do 1017 i=1,nn
      dss1(i,1)=ppxax(i)*pss1(i,1)+wxpx(i)*pss2(i,1)+c(i)
      dss1(i,2)=ppyay(i)*pss1(i,1)+wypy(i)*pss2(i,1)
      dss1(i,3)=ppzaz(i)*pss1(i,1)+wzpz(i)*pss2(i,1)
      dss1(i,4)=ppyay(i)*pss1(i,2)+wypy(i)*pss2(i,2)+c(i)
      dss1(i,5)=ppzaz(i)*pss1(i,2)+wzpz(i)*pss2(i,2)
 1017 dss1(i,6)=ppzaz(i)*pss1(i,3)+wzpz(i)*pss2(i,3)+c(i)

      do 1018 i=1,nn
 1018 c(i)=(sss2(i)-zetac1(i)*sss3(i))/(two*zzeta(i))

      do 1019 i=1,nn
      dss2(i,1)=ppxax(i)*pss2(i,1)+wxpx(i)*pss3(i,1)+c(i)
      dss2(i,2)=ppyay(i)*pss2(i,1)+wypy(i)*pss3(i,1)
      dss2(i,3)=ppzaz(i)*pss2(i,1)+wzpz(i)*pss3(i,1)
      dss2(i,4)=ppyay(i)*pss2(i,2)+wypy(i)*pss3(i,2)+c(i)
      dss2(i,5)=ppzaz(i)*pss2(i,2)+wzpz(i)*pss3(i,2)
 1019 dss2(i,6)=ppzaz(i)*pss2(i,3)+wzpz(i)*pss3(i,3)+c(i)

      do 1020 i=1,nn
 1020 c(i)=(sss3(i)-zetac1(i)*sss4(i))/(two*zzeta(i))

      do 1021 i=1,nn
      dss3(i,1)=ppxax(i)*pss3(i,1)+wxpx(i)*pss4(i,1)+c(i)
      dss3(i,2)=ppyay(i)*pss3(i,1)+wypy(i)*pss4(i,1)
      dss3(i,3)=ppzaz(i)*pss3(i,1)+wzpz(i)*pss4(i,1)
      dss3(i,4)=ppyay(i)*pss3(i,2)+wypy(i)*pss4(i,2)+c(i)
      dss3(i,5)=ppzaz(i)*pss3(i,2)+wzpz(i)*pss4(i,2)
 1021 dss3(i,6)=ppzaz(i)*pss3(i,3)+wzpz(i)*pss4(i,3)+c(i)

      do 1022 j=1,3
      do 1022 i=1,nn
 1022 c1(i,j)=(pss0(i,j)-zetac1(i)*pss1(i,j))/(two*zzeta(i))

      do 1023 i=1,nn
      fss0(i,01)=ppxax(i)*dss0(i,1)+wxpx(i)*dss1(i,1)+c1(i,1)+c1(i,1)
      fss0(i,02)=ppyay(i)*dss0(i,1)+wypy(i)*dss1(i,1)
      fss0(i,03)=ppzaz(i)*dss0(i,1)+wzpz(i)*dss1(i,1)
      fss0(i,04)=ppyay(i)*dss0(i,2)+wypy(i)*dss1(i,2)+c1(i,1)
      fss0(i,05)=ppzaz(i)*dss0(i,2)+wzpz(i)*dss1(i,2)
      fss0(i,06)=ppzaz(i)*dss0(i,3)+wzpz(i)*dss1(i,3)+c1(i,1)
      fss0(i,07)=ppyay(i)*dss0(i,4)+wypy(i)*dss1(i,4)+c1(i,2)+c1(i,2)
      fss0(i,08)=ppzaz(i)*dss0(i,4)+wzpz(i)*dss1(i,4)
      fss0(i,09)=ppzaz(i)*dss0(i,5)+wzpz(i)*dss1(i,5)+c1(i,2)
 1023 fss0(i,10)=ppzaz(i)*dss0(i,6)+wzpz(i)*dss1(i,6)+c1(i,3)+c1(i,3)

      do 1024 j=1,3
      do 1024 i=1,nn
 1024 c1(i,j)=(pss1(i,j)-zetac1(i)*pss2(i,j))/(two*zzeta(i))

      do 1025 i=1,nn
      fss1(i,01)=ppxax(i)*dss1(i,1)+wxpx(i)*dss2(i,1)+c1(i,1)+c1(i,1)
      fss1(i,02)=ppyay(i)*dss1(i,1)+wypy(i)*dss2(i,1)
      fss1(i,03)=ppzaz(i)*dss1(i,1)+wzpz(i)*dss2(i,1)
      fss1(i,04)=ppyay(i)*dss1(i,2)+wypy(i)*dss2(i,2)+c1(i,1)
      fss1(i,05)=ppzaz(i)*dss1(i,2)+wzpz(i)*dss2(i,2)
      fss1(i,06)=ppzaz(i)*dss1(i,3)+wzpz(i)*dss2(i,3)+c1(i,1)
      fss1(i,07)=ppyay(i)*dss1(i,4)+wypy(i)*dss2(i,4)+c1(i,2)+c1(i,2)
      fss1(i,08)=ppzaz(i)*dss1(i,4)+wzpz(i)*dss2(i,4)
      fss1(i,09)=ppzaz(i)*dss1(i,5)+wzpz(i)*dss2(i,5)+c1(i,2)
 1025 fss1(i,10)=ppzaz(i)*dss1(i,6)+wzpz(i)*dss2(i,6)+c1(i,3)+c1(i,3)

      do 1026 j=1,3
      do 1026 i=1,nn
 1026 c1(i,j)=(pss2(i,j)-zetac1(i)*pss3(i,j))/(two*zzeta(i))

      do 1027 i=1,nn
      fss2(i,01)=ppxax(i)*dss2(i,1)+wxpx(i)*dss3(i,1)+c1(i,1)+c1(i,1)
      fss2(i,02)=ppyay(i)*dss2(i,1)+wypy(i)*dss3(i,1)
      fss2(i,03)=ppzaz(i)*dss2(i,1)+wzpz(i)*dss3(i,1)
      fss2(i,04)=ppyay(i)*dss2(i,2)+wypy(i)*dss3(i,2)+c1(i,1)
      fss2(i,05)=ppzaz(i)*dss2(i,2)+wzpz(i)*dss3(i,2)
      fss2(i,06)=ppzaz(i)*dss2(i,3)+wzpz(i)*dss3(i,3)+c1(i,1)
      fss2(i,07)=ppyay(i)*dss2(i,4)+wypy(i)*dss3(i,4)+c1(i,2)+c1(i,2)
      fss2(i,08)=ppzaz(i)*dss2(i,4)+wzpz(i)*dss3(i,4)
      fss2(i,09)=ppzaz(i)*dss2(i,5)+wzpz(i)*dss3(i,5)+c1(i,2)
 1027 fss2(i,10)=ppzaz(i)*dss2(i,6)+wzpz(i)*dss3(i,6)+c1(i,3)+c1(i,3)

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

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

      do 1030 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)
 1030 psp0(i,9)=wzcz(i)*pss1(i,3)+c(i)

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

      do 1032 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)
 1032 psp1(i,9)=wzcz(i)*pss2(i,3)+c(i)

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

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

      do 1035 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)
 1035 psd0(i,18)=wzcz(i)*psp1(i,9)+c1(i,3)+c2(i,3)

      do 1036 j=1,3
      do 1036 i=1,nn
 1036 c1(i,j)=half*pss1(i,j)*zetarec(i)

      do 1037 i=1,nn
      dsp0(i,01)=wxcx(i)*dss1(i,1)+c1(i,1)*two
      dsp0(i,02)=wycy(i)*dss1(i,1)
      dsp0(i,03)=wzcz(i)*dss1(i,1)
      dsp0(i,04)=wxcx(i)*dss1(i,2)+c1(i,2)
      dsp0(i,05)=wycy(i)*dss1(i,2)+c1(i,1)
      dsp0(i,06)=wzcz(i)*dss1(i,2)
      dsp0(i,07)=wxcx(i)*dss1(i,3)+c1(i,3)
      dsp0(i,08)=wycy(i)*dss1(i,3)
      dsp0(i,09)=wzcz(i)*dss1(i,3)+c1(i,1)
      dsp0(i,10)=wxcx(i)*dss1(i,4)
      dsp0(i,11)=wycy(i)*dss1(i,4)+c1(i,2)*two
      dsp0(i,12)=wzcz(i)*dss1(i,4)
      dsp0(i,13)=wxcx(i)*dss1(i,5)
      dsp0(i,14)=wycy(i)*dss1(i,5)+c1(i,3)
      dsp0(i,15)=wzcz(i)*dss1(i,5)+c1(i,2)
      dsp0(i,16)=wxcx(i)*dss1(i,6)
      dsp0(i,17)=wycy(i)*dss1(i,6)
 1037 dsp0(i,18)=wzcz(i)*dss1(i,6)+c1(i,3)*two

      do 1038 j=1,3
      do 1038 i=1,nn
 1038 c1(i,j)=half*pss2(i,j)*zetarec(i)

      do 1039 i=1,nn
      dsp1(i,01)=wxcx(i)*dss2(i,1)+c1(i,1)*two
      dsp1(i,02)=wycy(i)*dss2(i,1)
      dsp1(i,03)=wzcz(i)*dss2(i,1)
      dsp1(i,04)=wxcx(i)*dss2(i,2)+c1(i,2)
      dsp1(i,05)=wycy(i)*dss2(i,2)+c1(i,1)
      dsp1(i,06)=wzcz(i)*dss2(i,2)
      dsp1(i,07)=wxcx(i)*dss2(i,3)+c1(i,3)
      dsp1(i,08)=wycy(i)*dss2(i,3)
      dsp1(i,09)=wzcz(i)*dss2(i,3)+c1(i,1)
      dsp1(i,10)=wxcx(i)*dss2(i,4)
      dsp1(i,11)=wycy(i)*dss2(i,4)+c1(i,2)*two
      dsp1(i,12)=wzcz(i)*dss2(i,4)
      dsp1(i,13)=wxcx(i)*dss2(i,5)
      dsp1(i,14)=wycy(i)*dss2(i,5)+c1(i,3)
      dsp1(i,15)=wzcz(i)*dss2(i,5)+c1(i,2)
      dsp1(i,16)=wxcx(i)*dss2(i,6)
      dsp1(i,17)=wycy(i)*dss2(i,6)
 1039 dsp1(i,18)=wzcz(i)*dss2(i,6)+c1(i,3)*two

      do 1040 j=1,6
      do 1040 i=1,nn
 1040 c1(i,j)=(dss0(i,j)-zetac2(i)*dss1(i,j))/(two*zetaux)

      do 1041 j=1,9
      do 1041 i=1,nn
 1041 c2(i,j)=half*psp1(i,j)*zetarec(i)

      do 1042 i=1,nn
      dsd0(i,01)=wxcx(i)*dsp1(i,01)+c1(i,1)+c2(i,1)*two
      dsd0(i,02)=wycy(i)*dsp1(i,01)
      dsd0(i,03)=wzcz(i)*dsp1(i,01)
      dsd0(i,04)=wycy(i)*dsp1(i,02)+c1(i,1)
      dsd0(i,05)=wzcz(i)*dsp1(i,02)
      dsd0(i,06)=wzcz(i)*dsp1(i,03)+c1(i,1)
      dsd0(i,07)=wxcx(i)*dsp1(i,04)+c1(i,2)+c2(i,4)
      dsd0(i,08)=wycy(i)*dsp1(i,04)        +c2(i,1)
      dsd0(i,09)=wzcz(i)*dsp1(i,04)
      dsd0(i,10)=wycy(i)*dsp1(i,05)+c1(i,2)+c2(i,2)
      dsd0(i,11)=wzcz(i)*dsp1(i,05)
      dsd0(i,12)=wzcz(i)*dsp1(i,06)+c1(i,2)
      dsd0(i,13)=wxcx(i)*dsp1(i,07)+c1(i,3)+c2(i,7)
      dsd0(i,14)=wycy(i)*dsp1(i,07)
      dsd0(i,15)=wzcz(i)*dsp1(i,07)        +c2(i,1)
      dsd0(i,16)=wycy(i)*dsp1(i,08)+c1(i,3)
      dsd0(i,17)=wzcz(i)*dsp1(i,08)        +c2(i,2)
      dsd0(i,18)=wzcz(i)*dsp1(i,09)+c1(i,3)+c2(i,3)
      dsd0(i,19)=wxcx(i)*dsp1(i,10)+c1(i,4)
      dsd0(i,20)=wycy(i)*dsp1(i,10)        +c2(i,4)*two
      dsd0(i,21)=wzcz(i)*dsp1(i,10)
      dsd0(i,22)=wycy(i)*dsp1(i,11)+c1(i,4)+c2(i,5)*two
      dsd0(i,23)=wzcz(i)*dsp1(i,11)
      dsd0(i,24)=wzcz(i)*dsp1(i,12)+c1(i,4)
      dsd0(i,25)=wxcx(i)*dsp1(i,13)+c1(i,5)
      dsd0(i,26)=wycy(i)*dsp1(i,13)        +c2(i,7)
      dsd0(i,27)=wzcz(i)*dsp1(i,13)        +c2(i,4)
      dsd0(i,28)=wycy(i)*dsp1(i,14)+c1(i,5)+c2(i,8)
      dsd0(i,29)=wzcz(i)*dsp1(i,14)        +c2(i,5)
      dsd0(i,30)=wzcz(i)*dsp1(i,15)+c1(i,5)+c2(i,6)
      dsd0(i,31)=wxcx(i)*dsp1(i,16)+c1(i,6)
      dsd0(i,32)=wycy(i)*dsp1(i,16)
      dsd0(i,33)=wzcz(i)*dsp1(i,16)        +c2(i,7)*two
      dsd0(i,34)=wycy(i)*dsp1(i,17)+c1(i,6)
      dsd0(i,35)=wzcz(i)*dsp1(i,17)        +c2(i,8)*two
 1042 dsd0(i,36)=wzcz(i)*dsp1(i,18)+c1(i,6)+c2(i,9)*two

      do 1043 j=1,6
      do 1043 i=1,nn
 1043 c1(i,j)=half*dss1(i,j)*zetarec(i)

      do 1044 i=1,nn
      fsp0(i,01)=wxcx(i)*fss1(i,01)+c1(i,1)*three
      fsp0(i,02)=wycy(i)*fss1(i,01)
      fsp0(i,03)=wzcz(i)*fss1(i,01)
      fsp0(i,04)=wxcx(i)*fss1(i,02)+c1(i,2)*two
      fsp0(i,05)=wycy(i)*fss1(i,02)+c1(i,1)
      fsp0(i,06)=wzcz(i)*fss1(i,02)
      fsp0(i,07)=wxcx(i)*fss1(i,03)+c1(i,3)*two
      fsp0(i,08)=wycy(i)*fss1(i,03)
      fsp0(i,09)=wzcz(i)*fss1(i,03)+c1(i,1)
      fsp0(i,10)=wxcx(i)*fss1(i,04)+c1(i,4)
      fsp0(i,11)=wycy(i)*fss1(i,04)+c1(i,2)*two
      fsp0(i,12)=wzcz(i)*fss1(i,04)
      fsp0(i,13)=wxcx(i)*fss1(i,05)+c1(i,5)
      fsp0(i,14)=wycy(i)*fss1(i,05)+c1(i,3)
      fsp0(i,15)=wzcz(i)*fss1(i,05)+c1(i,2)
      fsp0(i,16)=wxcx(i)*fss1(i,06)+c1(i,6)
      fsp0(i,17)=wycy(i)*fss1(i,06)
      fsp0(i,18)=wzcz(i)*fss1(i,06)+c1(i,3)*two
      fsp0(i,19)=wxcx(i)*fss1(i,07)
      fsp0(i,20)=wycy(i)*fss1(i,07)+c1(i,4)*three
      fsp0(i,21)=wzcz(i)*fss1(i,07)
      fsp0(i,22)=wxcx(i)*fss1(i,08)
      fsp0(i,23)=wycy(i)*fss1(i,08)+c1(i,5)*two
      fsp0(i,24)=wzcz(i)*fss1(i,08)+c1(i,4)
      fsp0(i,25)=wxcx(i)*fss1(i,09)
      fsp0(i,26)=wycy(i)*fss1(i,09)+c1(i,6)
      fsp0(i,27)=wzcz(i)*fss1(i,09)+c1(i,5)*two
      fsp0(i,28)=wxcx(i)*fss1(i,10)
      fsp0(i,29)=wycy(i)*fss1(i,10)
 1044 fsp0(i,30)=wzcz(i)*fss1(i,10)+c1(i,6)*three

      do 1045 j=1,6
      do 1045 i=1,nn
 1045 c1(i,j)=half*dss2(i,j)*zetarec(i)

      do 1046 i=1,nn
      fsp1(i,01)=wxcx(i)*fss2(i,01)+c1(i,1)*three
      fsp1(i,02)=wycy(i)*fss2(i,01)
      fsp1(i,03)=wzcz(i)*fss2(i,01)
      fsp1(i,04)=wxcx(i)*fss2(i,02)+c1(i,2)*two
      fsp1(i,05)=wycy(i)*fss2(i,02)+c1(i,1)
      fsp1(i,06)=wzcz(i)*fss2(i,02)
      fsp1(i,07)=wxcx(i)*fss2(i,03)+c1(i,3)*two
      fsp1(i,08)=wycy(i)*fss2(i,03)
      fsp1(i,09)=wzcz(i)*fss2(i,03)+c1(i,1)
      fsp1(i,10)=wxcx(i)*fss2(i,04)+c1(i,4)
      fsp1(i,11)=wycy(i)*fss2(i,04)+c1(i,2)*two
      fsp1(i,12)=wzcz(i)*fss2(i,04)
      fsp1(i,13)=wxcx(i)*fss2(i,05)+c1(i,5)
      fsp1(i,14)=wycy(i)*fss2(i,05)+c1(i,3)
      fsp1(i,15)=wzcz(i)*fss2(i,05)+c1(i,2)
      fsp1(i,16)=wxcx(i)*fss2(i,06)+c1(i,6)
      fsp1(i,17)=wycy(i)*fss2(i,06)
      fsp1(i,18)=wzcz(i)*fss2(i,06)+c1(i,3)*two
      fsp1(i,19)=wxcx(i)*fss2(i,07)
      fsp1(i,20)=wycy(i)*fss2(i,07)+c1(i,4)*three
      fsp1(i,21)=wzcz(i)*fss2(i,07)
      fsp1(i,22)=wxcx(i)*fss2(i,08)
      fsp1(i,23)=wycy(i)*fss2(i,08)+c1(i,5)*two
      fsp1(i,24)=wzcz(i)*fss2(i,08)+c1(i,4)
      fsp1(i,25)=wxcx(i)*fss2(i,09)
      fsp1(i,26)=wycy(i)*fss2(i,09)+c1(i,6)
      fsp1(i,27)=wzcz(i)*fss2(i,09)+c1(i,5)*two
      fsp1(i,28)=wxcx(i)*fss2(i,10)
      fsp1(i,29)=wycy(i)*fss2(i,10)
 1046 fsp1(i,30)=wzcz(i)*fss2(i,10)+c1(i,6)*three

      do 1047 j=1,10
      do 1047 i=1,nn
 1047 c1(i,j)=(fss0(i,j)-zetac2(i)*fss1(i,j))/(two*zetaux)

      do 1048 j=1,18
      do 1048 i=1,nn
 1048 c2(i,j)=half*dsp1(i,j)*zetarec(i)

      do 1049 i=1,nn
      fsd0(i,01)=wxcx(i)*fsp1(i,01)+c1(i,01)+c2(i,01)*three
      fsd0(i,02)=wycy(i)*fsp1(i,01)
      fsd0(i,03)=wzcz(i)*fsp1(i,01)
      fsd0(i,04)=wycy(i)*fsp1(i,02)+c1(i,01)
      fsd0(i,05)=wzcz(i)*fsp1(i,02)
      fsd0(i,06)=wzcz(i)*fsp1(i,03)+c1(i,01)
      fsd0(i,07)=wxcx(i)*fsp1(i,04)+c1(i,02)+c2(i,04)*two
      fsd0(i,08)=wycy(i)*fsp1(i,04)         +c2(i,01)
      fsd0(i,09)=wzcz(i)*fsp1(i,04)
      fsd0(i,10)=wycy(i)*fsp1(i,05)+c1(i,02)+c2(i,02)
      fsd0(i,11)=wzcz(i)*fsp1(i,05)
      fsd0(i,12)=wzcz(i)*fsp1(i,06)+c1(i,02)
      fsd0(i,13)=wxcx(i)*fsp1(i,07)+c1(i,03)+c2(i,07)*two
      fsd0(i,14)=wycy(i)*fsp1(i,07)
      fsd0(i,15)=wzcz(i)*fsp1(i,07)         +c2(i,01)
      fsd0(i,16)=wycy(i)*fsp1(i,08)+c1(i,03)
      fsd0(i,17)=wzcz(i)*fsp1(i,08)         +c2(i,02)
      fsd0(i,18)=wzcz(i)*fsp1(i,09)+c1(i,03)+c2(i,03)
      fsd0(i,19)=wxcx(i)*fsp1(i,10)+c1(i,04)+c2(i,10)
      fsd0(i,20)=wycy(i)*fsp1(i,10)         +c2(i,04)*two
      fsd0(i,21)=wzcz(i)*fsp1(i,10)
      fsd0(i,22)=wycy(i)*fsp1(i,11)+c1(i,04)+c2(i,05)*two
      fsd0(i,23)=wzcz(i)*fsp1(i,11)
      fsd0(i,24)=wzcz(i)*fsp1(i,12)+c1(i,04)
      fsd0(i,25)=wxcx(i)*fsp1(i,13)+c1(i,05)+c2(i,13)
      fsd0(i,26)=wycy(i)*fsp1(i,13)         +c2(i,07)
      fsd0(i,27)=wzcz(i)*fsp1(i,13)         +c2(i,04)
      fsd0(i,28)=wycy(i)*fsp1(i,14)+c1(i,05)+c2(i,08)
      fsd0(i,29)=wzcz(i)*fsp1(i,14)         +c2(i,05)
      fsd0(i,30)=wzcz(i)*fsp1(i,15)+c1(i,05)+c2(i,06)
      fsd0(i,31)=wxcx(i)*fsp1(i,16)+c1(i,06)+c2(i,16)
      fsd0(i,32)=wycy(i)*fsp1(i,16)
      fsd0(i,33)=wzcz(i)*fsp1(i,16)         +c2(i,07)*two
      fsd0(i,34)=wycy(i)*fsp1(i,17)+c1(i,06)
      fsd0(i,35)=wzcz(i)*fsp1(i,17)         +c2(i,08)*two
      fsd0(i,36)=wzcz(i)*fsp1(i,18)+c1(i,06)+c2(i,09)*two
      fsd0(i,37)=wxcx(i)*fsp1(i,19)+c1(i,07)
      fsd0(i,38)=wycy(i)*fsp1(i,19)         +c2(i,10)*three
      fsd0(i,39)=wzcz(i)*fsp1(i,19)
      fsd0(i,40)=wycy(i)*fsp1(i,20)+c1(i,07)+c2(i,11)*three
      fsd0(i,41)=wzcz(i)*fsp1(i,20)
      fsd0(i,42)=wzcz(i)*fsp1(i,21)+c1(i,07)
      fsd0(i,43)=wxcx(i)*fsp1(i,22)+c1(i,08)
      fsd0(i,44)=wycy(i)*fsp1(i,22)         +c2(i,13)*two
      fsd0(i,45)=wzcz(i)*fsp1(i,22)         +c2(i,10)
      fsd0(i,46)=wycy(i)*fsp1(i,23)+c1(i,08)+c2(i,14)*two
      fsd0(i,47)=wzcz(i)*fsp1(i,23)         +c2(i,11)
      fsd0(i,48)=wzcz(i)*fsp1(i,24)+c1(i,08)+c2(i,12)
      fsd0(i,49)=wxcx(i)*fsp1(i,25)+c1(i,09)
      fsd0(i,50)=wycy(i)*fsp1(i,25)         +c2(i,16)
      fsd0(i,51)=wzcz(i)*fsp1(i,25)         +c2(i,13)*two
      fsd0(i,52)=wycy(i)*fsp1(i,26)+c1(i,09)+c2(i,17)
      fsd0(i,53)=wzcz(i)*fsp1(i,26)         +c2(i,14)*two
      fsd0(i,54)=wzcz(i)*fsp1(i,27)+c1(i,09)+c2(i,15)*two
      fsd0(i,55)=wxcx(i)*fsp1(i,28)+c1(i,10)
      fsd0(i,56)=wycy(i)*fsp1(i,28)
      fsd0(i,57)=wzcz(i)*fsp1(i,28)         +c2(i,16)*three
      fsd0(i,58)=wycy(i)*fsp1(i,29)+c1(i,10)
      fsd0(i,59)=wzcz(i)*fsp1(i,29)         +c2(i,17)*three
 1049 fsd0(i,60)=wzcz(i)*fsp1(i,30)+c1(i,10)+c2(i,18)*three

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

      do 1051 j=1,3
      do 1051 i=1,nn
 1051 psscond(iicontr(i),j)=psscond(iicontr(i),j)
     &                     -pss0(i,j)*cdfits

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

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

      do 1054 i=1,m
      grad(i,1)=psscond(i,1)*density(inddyy(i,1))*two
      grad(i,2)=psscond(i,2)*density(inddyy(i,2))
      grad(i,3)=psscond(i,3)*density(inddyy(i,3))
      grad(i,4)=psscond(i,1)*density(inddyy(i,2))
      grad(i,5)=psscond(i,2)*density(inddyy(i,4))*two
      grad(i,6)=psscond(i,3)*density(inddyy(i,5))
      grad(i,7)=psscond(i,1)*density(inddyy(i,3))
      grad(i,8)=psscond(i,2)*density(inddyy(i,5))
 1054 grad(i,9)=psscond(i,3)*density(inddyy(i,6))*two

      do 1055 k=1,3
      do 1055 j=1,3
      do 1055 i=1,m
 1055 gradient(k,nc)=gradient(k,nc)-grad(i,3*(k-1)+j)

      do 1056 k=1,3
      do 1056 j=1,3
      do 1056 i=1,m
 1056 gradient(k,icentr(i,1))=gradient(k,icentr(i,1))+grad(i,3*(k-1)+j)

      do 1057 j=1,10
      do 1057 i=1,m
 1057 fsscond(i,j)=zero

      do 1058 j=1,10
      do 1058 i=1,nn
 1058 fsscond(iicontr(i),j)=fsscond(iicontr(i),j)
     &                     +fss0(i,j)*ttwozetaa(i)*cdfits

      do 1059 k=1,3
      do 1059 j=1,10
      do 1059 i=1,nn
 1059 fsscond(iicontr(i),j)=fsscond(iicontr(i),j)
     &                     +fsp0(i,3*(j-1)+k)*ttwozetaa(i)*cdfitp(k)

      do 1060 k=1,6
      do 1060 j=1,10
      do 1060 i=1,nn
 1060 fsscond(iicontr(i),j)=fsscond(iicontr(i),j)
     &                     +fsd0(i,6*(j-1)+k)*ttwozetaa(i)*cdfitd(k)

      do 1061 i=1,m
      grad(i,01)=fsscond(i,01)*density(inddyy(i,1))
      grad(i,02)=fsscond(i,02)*density(inddyy(i,2))
      grad(i,03)=fsscond(i,03)*density(inddyy(i,3))
      grad(i,04)=fsscond(i,04)*density(inddyy(i,4))
      grad(i,05)=fsscond(i,05)*density(inddyy(i,5))
      grad(i,06)=fsscond(i,06)*density(inddyy(i,6))
      grad(i,07)=fsscond(i,02)*density(inddyy(i,1))
      grad(i,08)=fsscond(i,04)*density(inddyy(i,2))
      grad(i,09)=fsscond(i,05)*density(inddyy(i,3))
      grad(i,10)=fsscond(i,07)*density(inddyy(i,4))
      grad(i,11)=fsscond(i,08)*density(inddyy(i,5))
      grad(i,12)=fsscond(i,09)*density(inddyy(i,6))
      grad(i,13)=fsscond(i,03)*density(inddyy(i,1))
      grad(i,14)=fsscond(i,05)*density(inddyy(i,2))
      grad(i,15)=fsscond(i,06)*density(inddyy(i,3))
      grad(i,16)=fsscond(i,08)*density(inddyy(i,4))
      grad(i,17)=fsscond(i,09)*density(inddyy(i,5))
 1061 grad(i,18)=fsscond(i,10)*density(inddyy(i,6))

      do 1062 k=1,3
      do 1062 j=1,6
      do 1062 i=1,m
 1062 gradient(k,nc)=gradient(k,nc)-grad(i,6*(k-1)+j)

      do 1063 k=1,3
      do 1063 j=1,6
      do 1063 i=1,m
 1063 gradient(k,icentr(i,1))=gradient(k,icentr(i,1))+grad(i,6*(k-1)+j)

      do 1064 j=1,6
      do 1064 i=1,m
 1064 dsscond(i,j)=zero

      do 1065 j=1,10
      do 1065 i=1,m
 1065 fsscond(i,j)=zero

      do 1066 j=1,6
      do 1066 i=1,nn
 1066 dsscond(iicontr(i),j)=dsscond(iicontr(i),j)
     &                     +dss0(i,j)*ttwozetab(i)*cdfits

      do 1067 k=1,3
      do 1067 j=1,6
      do 1067 i=1,nn
 1067 dsscond(iicontr(i),j)=dsscond(iicontr(i),j)
     &                     +dsp0(i,3*(j-1)+k)*ttwozetab(i)*cdfitp(k)

      do 1068 k=1,6
      do 1068 j=1,6
      do 1068 i=1,nn
 1068 dsscond(iicontr(i),j)=dsscond(iicontr(i),j)
     &                     +dsd0(i,6*(j-1)+k)*ttwozetab(i)*cdfitd(k)

      do 1069 j=1,10
      do 1069 i=1,nn
 1069 fsscond(iicontr(i),j)=fsscond(iicontr(i),j)
     &                     +fss0(i,j)*ttwozetab(i)*cdfits

      do 1070 k=1,3
      do 1070 j=1,10
      do 1070 i=1,nn
 1070 fsscond(iicontr(i),j)=fsscond(iicontr(i),j)
     &                     +fsp0(i,3*(j-1)+k)*ttwozetab(i)*cdfitp(k)

      do 1071 k=1,6
      do 1071 j=1,10
      do 1071 i=1,nn
 1071 fsscond(iicontr(i),j)=fsscond(iicontr(i),j)
     &                     +fsd0(i,6*(j-1)+k)*ttwozetab(i)*cdfitd(k)

      do 1072 i=1,m
      dpscond(i,01)=axbx(i)*dsscond(i,1)+fsscond(i,01)
      dpscond(i,02)=ayby(i)*dsscond(i,1)+fsscond(i,02)
      dpscond(i,03)=azbz(i)*dsscond(i,1)+fsscond(i,03)
      dpscond(i,04)=axbx(i)*dsscond(i,2)+fsscond(i,02)
      dpscond(i,05)=ayby(i)*dsscond(i,2)+fsscond(i,04)
      dpscond(i,06)=azbz(i)*dsscond(i,2)+fsscond(i,05)
      dpscond(i,07)=axbx(i)*dsscond(i,3)+fsscond(i,03)
      dpscond(i,08)=ayby(i)*dsscond(i,3)+fsscond(i,05)
      dpscond(i,09)=azbz(i)*dsscond(i,3)+fsscond(i,06)
      dpscond(i,10)=axbx(i)*dsscond(i,4)+fsscond(i,04)
      dpscond(i,11)=ayby(i)*dsscond(i,4)+fsscond(i,07)
      dpscond(i,12)=azbz(i)*dsscond(i,4)+fsscond(i,08)
      dpscond(i,13)=axbx(i)*dsscond(i,5)+fsscond(i,05)
      dpscond(i,14)=ayby(i)*dsscond(i,5)+fsscond(i,08)
      dpscond(i,15)=azbz(i)*dsscond(i,5)+fsscond(i,09)
      dpscond(i,16)=axbx(i)*dsscond(i,6)+fsscond(i,06)
      dpscond(i,17)=ayby(i)*dsscond(i,6)+fsscond(i,09)
 1072 dpscond(i,18)=azbz(i)*dsscond(i,6)+fsscond(i,10)

      do 1073 i=1,m
      grad(i,01)=dpscond(i,01)*density(inddyy(i,1))
      grad(i,02)=dpscond(i,04)*density(inddyy(i,2))
      grad(i,03)=dpscond(i,07)*density(inddyy(i,3))
      grad(i,04)=dpscond(i,10)*density(inddyy(i,4))
      grad(i,05)=dpscond(i,13)*density(inddyy(i,5))
      grad(i,06)=dpscond(i,16)*density(inddyy(i,6))
      grad(i,07)=dpscond(i,02)*density(inddyy(i,1))
      grad(i,08)=dpscond(i,05)*density(inddyy(i,2))
      grad(i,09)=dpscond(i,08)*density(inddyy(i,3))
      grad(i,10)=dpscond(i,11)*density(inddyy(i,4))
      grad(i,11)=dpscond(i,14)*density(inddyy(i,5))
      grad(i,12)=dpscond(i,17)*density(inddyy(i,6))
      grad(i,13)=dpscond(i,03)*density(inddyy(i,1))
      grad(i,14)=dpscond(i,06)*density(inddyy(i,2))
      grad(i,15)=dpscond(i,09)*density(inddyy(i,3))
      grad(i,16)=dpscond(i,12)*density(inddyy(i,4))
      grad(i,17)=dpscond(i,15)*density(inddyy(i,5))
 1073 grad(i,18)=dpscond(i,18)*density(inddyy(i,6))

      do 1074 k=1,3
      do 1074 j=1,6
      do 1074 i=1,m
 1074 gradient(k,nc)=gradient(k,nc)-grad(i,6*(k-1)+j)

      do 1075 k=1,3
      do 1075 j=1,6
      do 1075 i=1,m
 1075 gradient(k,icentr(i,2))=gradient(k,icentr(i,2))+grad(i,6*(k-1)+j)

 1001 continue

      return
      end
