      subroutine gradxxdp3(m,n,nn,nc,inddyy,icentr,iicontr,cdfits,
     &                     cdfitp,cdfitd,gradient,density,ttwozetab,
     &                     dss0,dsp0,dsd0,fss0,fsp0,fsd0,gss0,gsp0,gsd0,
     &                     axbx,ayby,azbz,dsscond,dspcond,dsdcond,
     &                     fsscond,fspcond,fsdcond,gsscond,gspcond,
     &                     gsdcond,dpscond,dppcond,dpdcond,fpscond,
     &                     fppcond,fpdcond,ddscond,ddpcond,dddcond,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 inddyy(m,*),icentr(m,*),iicontr(*)

      dimension cdfitp(*),cdfitd(*),gradient(3,*),density(*),
     &          ttwozetab(*),dss0(n,*),dsp0(n,*),dsd0(n,*),fss0(n,*),
     &          fsp0(n,*),fsd0(n,*),gss0(n,*),gsp0(n,*),gsd0(n,*),
     &          axbx(*),ayby(*),azbz(*),dsscond(m,*),dspcond(m,*),
     &          dsdcond(m,*),fsscond(m,*),fspcond(m,*),fsdcond(m,*),
     &          gsscond(m,*),gspcond(m,*),gsdcond(m,*),dpscond(m,*),
     &          dppcond(m,*),dpdcond(m,*),fpscond(m,*),fppcond(m,*),
     &          fpdcond(m,*),ddscond(m,*),ddpcond(m,*),dddcond(m,*),
     &          grad(m,*)

      data zero/0.0d0/

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

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

      do 1003 j=1,15
      do 1003 i=1,m
 1003 gsscond(i,j)=zero

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

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

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

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

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

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

      do 1010 j=1,15
      do 1010 i=1,nn
 1010 gsscond(iicontr(i),j)=gsscond(iicontr(i),j)
     &                     +gss0(i,j)*ttwozetab(i)*cdfits

      do 1011 k=1,3
      do 1011 j=1,15
      do 1011 i=1,nn
 1011 gsscond(iicontr(i),j)=gsscond(iicontr(i),j)
     &                     +gsp0(i,3*(j-1)+k)*ttwozetab(i)*cdfitp(k)

      do 1012 k=1,6
      do 1012 j=1,15
      do 1012 i=1,nn
 1012 gsscond(iicontr(i),j)=gsscond(iicontr(i),j)
     &                     +gsd0(i,6*(j-1)+k)*ttwozetab(i)*cdfitd(k)

      do 1013 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)
 1013 dpscond(i,18)=azbz(i)*dsscond(i,6)+fsscond(i,10)

      do 1014 i=1,m
      fpscond(i,01)=axbx(i)*fsscond(i,01)+gsscond(i,01)
      fpscond(i,02)=ayby(i)*fsscond(i,01)+gsscond(i,02)
      fpscond(i,03)=azbz(i)*fsscond(i,01)+gsscond(i,03)
      fpscond(i,04)=axbx(i)*fsscond(i,02)+gsscond(i,02)
      fpscond(i,05)=ayby(i)*fsscond(i,02)+gsscond(i,04)
      fpscond(i,06)=azbz(i)*fsscond(i,02)+gsscond(i,05)
      fpscond(i,07)=axbx(i)*fsscond(i,03)+gsscond(i,03)
      fpscond(i,08)=ayby(i)*fsscond(i,03)+gsscond(i,05)
      fpscond(i,09)=azbz(i)*fsscond(i,03)+gsscond(i,06)
      fpscond(i,10)=axbx(i)*fsscond(i,04)+gsscond(i,04)
      fpscond(i,11)=ayby(i)*fsscond(i,04)+gsscond(i,07)
      fpscond(i,12)=azbz(i)*fsscond(i,04)+gsscond(i,08)
      fpscond(i,13)=axbx(i)*fsscond(i,05)+gsscond(i,05)
      fpscond(i,14)=ayby(i)*fsscond(i,05)+gsscond(i,08)
      fpscond(i,15)=azbz(i)*fsscond(i,05)+gsscond(i,09)
      fpscond(i,16)=axbx(i)*fsscond(i,06)+gsscond(i,06)
      fpscond(i,17)=ayby(i)*fsscond(i,06)+gsscond(i,09)
      fpscond(i,18)=azbz(i)*fsscond(i,06)+gsscond(i,10)
      fpscond(i,19)=axbx(i)*fsscond(i,07)+gsscond(i,07)
      fpscond(i,20)=ayby(i)*fsscond(i,07)+gsscond(i,11)
      fpscond(i,21)=azbz(i)*fsscond(i,07)+gsscond(i,12)
      fpscond(i,22)=axbx(i)*fsscond(i,08)+gsscond(i,08)
      fpscond(i,23)=ayby(i)*fsscond(i,08)+gsscond(i,12)
      fpscond(i,24)=azbz(i)*fsscond(i,08)+gsscond(i,13)
      fpscond(i,25)=axbx(i)*fsscond(i,09)+gsscond(i,09)
      fpscond(i,26)=ayby(i)*fsscond(i,09)+gsscond(i,13)
      fpscond(i,27)=azbz(i)*fsscond(i,09)+gsscond(i,14)
      fpscond(i,28)=axbx(i)*fsscond(i,10)+gsscond(i,10)
      fpscond(i,29)=ayby(i)*fsscond(i,10)+gsscond(i,14)
 1014 fpscond(i,30)=azbz(i)*fsscond(i,10)+gsscond(i,15)

      do 1015 i=1,m
      ddscond(i,01)=axbx(i)*dpscond(i,01)+fpscond(i,01)
      ddscond(i,02)=ayby(i)*dpscond(i,01)+fpscond(i,04)
      ddscond(i,03)=azbz(i)*dpscond(i,01)+fpscond(i,07)
      ddscond(i,04)=ayby(i)*dpscond(i,02)+fpscond(i,05)
      ddscond(i,05)=azbz(i)*dpscond(i,02)+fpscond(i,08)
      ddscond(i,06)=azbz(i)*dpscond(i,03)+fpscond(i,09)
      ddscond(i,07)=axbx(i)*dpscond(i,04)+fpscond(i,04)
      ddscond(i,08)=ayby(i)*dpscond(i,04)+fpscond(i,10)
      ddscond(i,09)=azbz(i)*dpscond(i,04)+fpscond(i,13)
      ddscond(i,10)=ayby(i)*dpscond(i,05)+fpscond(i,11)
      ddscond(i,11)=azbz(i)*dpscond(i,05)+fpscond(i,14)
      ddscond(i,12)=azbz(i)*dpscond(i,06)+fpscond(i,15)
      ddscond(i,13)=axbx(i)*dpscond(i,07)+fpscond(i,07)
      ddscond(i,14)=ayby(i)*dpscond(i,07)+fpscond(i,13)
      ddscond(i,15)=azbz(i)*dpscond(i,07)+fpscond(i,16)
      ddscond(i,16)=ayby(i)*dpscond(i,08)+fpscond(i,14)
      ddscond(i,17)=azbz(i)*dpscond(i,08)+fpscond(i,17)
      ddscond(i,18)=azbz(i)*dpscond(i,09)+fpscond(i,18)
      ddscond(i,19)=axbx(i)*dpscond(i,10)+fpscond(i,10)
      ddscond(i,20)=ayby(i)*dpscond(i,10)+fpscond(i,19)
      ddscond(i,21)=azbz(i)*dpscond(i,10)+fpscond(i,22)
      ddscond(i,22)=ayby(i)*dpscond(i,11)+fpscond(i,20)
      ddscond(i,23)=azbz(i)*dpscond(i,11)+fpscond(i,23)
      ddscond(i,24)=azbz(i)*dpscond(i,12)+fpscond(i,24)
      ddscond(i,25)=axbx(i)*dpscond(i,13)+fpscond(i,13)
      ddscond(i,26)=ayby(i)*dpscond(i,13)+fpscond(i,22)
      ddscond(i,27)=azbz(i)*dpscond(i,13)+fpscond(i,25)
      ddscond(i,28)=ayby(i)*dpscond(i,14)+fpscond(i,23)
      ddscond(i,29)=azbz(i)*dpscond(i,14)+fpscond(i,26)
      ddscond(i,30)=azbz(i)*dpscond(i,15)+fpscond(i,27)
      ddscond(i,31)=axbx(i)*dpscond(i,16)+fpscond(i,16)
      ddscond(i,32)=ayby(i)*dpscond(i,16)+fpscond(i,25)
      ddscond(i,33)=azbz(i)*dpscond(i,16)+fpscond(i,28)
      ddscond(i,34)=ayby(i)*dpscond(i,17)+fpscond(i,26)
      ddscond(i,35)=azbz(i)*dpscond(i,17)+fpscond(i,29)
 1015 ddscond(i,36)=azbz(i)*dpscond(i,18)+fpscond(i,30)

      do 1016 i=1,m
      grad(i,01)=ddscond(i,01)*density(inddyy(i,01))
      grad(i,02)=ddscond(i,02)*density(inddyy(i,02))
      grad(i,03)=ddscond(i,03)*density(inddyy(i,03))
      grad(i,04)=ddscond(i,07)*density(inddyy(i,04))
      grad(i,05)=ddscond(i,08)*density(inddyy(i,05))
      grad(i,06)=ddscond(i,09)*density(inddyy(i,06))
      grad(i,07)=ddscond(i,13)*density(inddyy(i,07))
      grad(i,08)=ddscond(i,14)*density(inddyy(i,08))
      grad(i,09)=ddscond(i,15)*density(inddyy(i,09))
      grad(i,10)=ddscond(i,19)*density(inddyy(i,10))
      grad(i,11)=ddscond(i,20)*density(inddyy(i,11))
      grad(i,12)=ddscond(i,21)*density(inddyy(i,12))
      grad(i,13)=ddscond(i,25)*density(inddyy(i,13))
      grad(i,14)=ddscond(i,26)*density(inddyy(i,14))
      grad(i,15)=ddscond(i,27)*density(inddyy(i,15))
      grad(i,16)=ddscond(i,31)*density(inddyy(i,16))
      grad(i,17)=ddscond(i,32)*density(inddyy(i,17))
      grad(i,18)=ddscond(i,33)*density(inddyy(i,18))
      grad(i,19)=ddscond(i,02)*density(inddyy(i,01))
      grad(i,20)=ddscond(i,04)*density(inddyy(i,02))
      grad(i,21)=ddscond(i,05)*density(inddyy(i,03))
      grad(i,22)=ddscond(i,08)*density(inddyy(i,04))
      grad(i,23)=ddscond(i,10)*density(inddyy(i,05))
      grad(i,24)=ddscond(i,11)*density(inddyy(i,06))
      grad(i,25)=ddscond(i,14)*density(inddyy(i,07))
      grad(i,26)=ddscond(i,16)*density(inddyy(i,08))
      grad(i,27)=ddscond(i,17)*density(inddyy(i,09))
      grad(i,28)=ddscond(i,20)*density(inddyy(i,10))
      grad(i,29)=ddscond(i,22)*density(inddyy(i,11))
      grad(i,30)=ddscond(i,23)*density(inddyy(i,12))
      grad(i,31)=ddscond(i,26)*density(inddyy(i,13))
      grad(i,32)=ddscond(i,28)*density(inddyy(i,14))
      grad(i,33)=ddscond(i,29)*density(inddyy(i,15))
      grad(i,34)=ddscond(i,32)*density(inddyy(i,16))
      grad(i,35)=ddscond(i,34)*density(inddyy(i,17))
      grad(i,36)=ddscond(i,35)*density(inddyy(i,18))
      grad(i,37)=ddscond(i,03)*density(inddyy(i,01))
      grad(i,38)=ddscond(i,05)*density(inddyy(i,02))
      grad(i,39)=ddscond(i,06)*density(inddyy(i,03))
      grad(i,40)=ddscond(i,09)*density(inddyy(i,04))
      grad(i,41)=ddscond(i,11)*density(inddyy(i,05))
      grad(i,42)=ddscond(i,12)*density(inddyy(i,06))
      grad(i,43)=ddscond(i,15)*density(inddyy(i,07))
      grad(i,44)=ddscond(i,17)*density(inddyy(i,08))
      grad(i,45)=ddscond(i,18)*density(inddyy(i,09))
      grad(i,46)=ddscond(i,21)*density(inddyy(i,10))
      grad(i,47)=ddscond(i,23)*density(inddyy(i,11))
      grad(i,48)=ddscond(i,24)*density(inddyy(i,12))
      grad(i,49)=ddscond(i,27)*density(inddyy(i,13))
      grad(i,50)=ddscond(i,29)*density(inddyy(i,14))
      grad(i,51)=ddscond(i,30)*density(inddyy(i,15))
      grad(i,52)=ddscond(i,33)*density(inddyy(i,16))
      grad(i,53)=ddscond(i,35)*density(inddyy(i,17))
 1016 grad(i,54)=ddscond(i,36)*density(inddyy(i,18))

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

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

      return
      end
