      subroutine gradxxdd2(m,n,nn,nc,iminj,inddyy,icentr,iicontr,cdfits,
     &                     cdfitp,cdfitd,gradient,density,pss0,psp0,
     &                     psd0,dss0,dsp0,dsd0,fss0,fsp0,fsd0,axbx,ayby,
     &                     azbz,psscond,pspcond,psdcond,dsscond,dspcond,
     &                     dsdcond,fsscond,fspcond,fsdcond,ppscond,
     &                     pppcond,ppdcond,dpscond,dppcond,dpdcond,
     &                     pdscond,pdpcond,pddcond,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 iminj(*),inddyy(m,*),icentr(m,*),iicontr(*)

      dimension cdfitp(*),cdfitd(*),gradient(3,*),density(*),pss0(n,*),
     &          psp0(n,*),psd0(n,*),dss0(n,*),dsp0(n,*),dsd0(n,*),
     &          fss0(n,*),fsp0(n,*),fsd0(n,*),axbx(*),ayby(*),azbz(*),
     &          psscond(m,*),pspcond(m,*),psdcond(m,*),dsscond(m,*),
     &          dspcond(m,*),dsdcond(m,*),fsscond(m,*),fspcond(m,*),
     &          fsdcond(m,*),ppscond(m,*),pppcond(m,*),ppdcond(m,*),
     &          dpscond(m,*),dppcond(m,*),dpdcond(m,*),pdscond(m,*),
     &          pdpcond(m,*),pddcond(m,*),grad(m,*)

      data zero,two/0.0d0,2.0d0/

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

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

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

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

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

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

      do 1007 j=1,6
      do 1007 i=1,nn
 1007 dsscond(iicontr(i),j)=dsscond(iicontr(i),j)
     &                     -dss0(i,j)*cdfits

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

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

      do 1010 j=1,10
      do 1010 i=1,nn
 1010 fsscond(iicontr(i),j)=fsscond(iicontr(i),j)
     &                     -fss0(i,j)*cdfits

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

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

      do 1013 i=1,m
      ppscond(i,1)=axbx(i)*psscond(i,1)+dsscond(i,1)
      ppscond(i,2)=ayby(i)*psscond(i,1)+dsscond(i,2)
      ppscond(i,3)=azbz(i)*psscond(i,1)+dsscond(i,3)
      ppscond(i,4)=axbx(i)*psscond(i,2)+dsscond(i,2)
      ppscond(i,5)=ayby(i)*psscond(i,2)+dsscond(i,4)
      ppscond(i,6)=azbz(i)*psscond(i,2)+dsscond(i,5)
      ppscond(i,7)=axbx(i)*psscond(i,3)+dsscond(i,3)
      ppscond(i,8)=ayby(i)*psscond(i,3)+dsscond(i,5)
 1013 ppscond(i,9)=azbz(i)*psscond(i,3)+dsscond(i,6)

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

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

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

      do 1017 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,02)=zero
                              grad(i,03)=zero
                              grad(i,04)=zero
                              grad(i,05)=zero
                              grad(i,06)=zero
                              grad(i,09)=zero
                              grad(i,10)=zero
                              grad(i,11)=zero
                              grad(i,12)=zero
                              grad(i,16)=zero
                              grad(i,17)=zero
                              grad(i,18)=zero
                              grad(i,21)=zero
                              grad(i,22)=zero
                              grad(i,23)=zero
                              grad(i,24)=zero
                              grad(i,29)=zero
                              grad(i,30)=zero
                              grad(i,36)=zero
                              grad(i,40)=zero
                              grad(i,41)=zero
                              grad(i,42)=zero
                              grad(i,48)=zero
                        endif
 1017 continue

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

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

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

      do 1021 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,02)=zero
                              grad(i,03)=zero
                              grad(i,06)=zero
                              grad(i,19)=zero
                              grad(i,20)=zero
                              grad(i,21)=zero
                              grad(i,23)=zero
                              grad(i,24)=zero
                              grad(i,26)=zero
                              grad(i,27)=zero
                              grad(i,30)=zero
                              grad(i,37)=zero
                              grad(i,38)=zero
                              grad(i,39)=zero
                              grad(i,40)=zero
                              grad(i,41)=zero
                              grad(i,42)=zero
                              grad(i,44)=zero
                              grad(i,45)=zero
                              grad(i,47)=zero
                              grad(i,48)=zero
                              grad(i,51)=zero
                        endif
 1021 continue

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

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

      return
      end
