      subroutine gradovdd(ndim,m,n,iminj,icontr,inddyy,icentr,density,
     &                    overlap,gradient,zeta,twozetaa,twozetab,eta,
     &                    eta2,base,etarr,pxax,pyay,pzaz,pxbx,pyby,pzbz,
     &                    twozeta,ss0,ss,ps0,ps,sp,ds0,ds,sd,fs0,fs,dp,
     &                    pd,fp,fd,c,c1,c2,pdcond,fdcond,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(*),icontr(*),inddyy(m,*),icentr(m,*)

      dimension density(*),overlap(ndim,*),gradient(3,*),zeta(*),
     &          twozetaa(*),twozetab(*),eta(*),eta2(*),base(*),etarr(*),
     &          pxax(*),pyay(*),pzaz(*),pxbx(*),pyby(*),pzbz(*),
     &          twozeta(*),ss0(*),ss(*),ps0(n,*),ps(n,*),sp(n,*),
     &          ds0(n,*),ds(n,*),sd(n,*),fs0(n,*),fs(n,*),dp(n,*),
     &          pd(n,*),fp(n,*),fd(n,*),c(*),c1(n,*),c2(n,*),
     &          pdcond(m,*),fdcond(m,*),grad(m,*)

      data zero,thrhlf,two,three,pi/0.0d0,1.5d0,
     &                              2.0d0,3.0d0,3.1415926535898d0/

      do 1001 i=1,n
 1001 twozeta(i)=two*zeta(i)

      do 1002 i=1,n
 1002 ss0(i)=((pi/zeta(i))**thrhlf)*base(i)

      do 1003 i=1,n
      ps0(i,1)=pxax(i)*ss0(i)
      ps0(i,2)=pyay(i)*ss0(i)
 1003 ps0(i,3)=pzaz(i)*ss0(i)

      do 1004 i=1,n
      sp(i,1)=pxbx(i)*ss0(i)
      sp(i,2)=pyby(i)*ss0(i)
 1004 sp(i,3)=pzbz(i)*ss0(i)

      do 1005 i=1,n
 1005 c(i)=ss0(i)/twozeta(i)

      do 1006 i=1,n
      ds0(i,1)=pxax(i)*ps0(i,1)+c(i)
      ds0(i,2)=pyay(i)*ps0(i,1)
      ds0(i,3)=pzaz(i)*ps0(i,1)
      ds0(i,4)=pyay(i)*ps0(i,2)+c(i)
      ds0(i,5)=pzaz(i)*ps0(i,2)
 1006 ds0(i,6)=pzaz(i)*ps0(i,3)+c(i)

      do 1007 i=1,n
      sd(i,1)=pxbx(i)*sp(i,1)+c(i)
      sd(i,2)=pyby(i)*sp(i,1)
      sd(i,3)=pzbz(i)*sp(i,1)
      sd(i,4)=pyby(i)*sp(i,2)+c(i)
      sd(i,5)=pzbz(i)*sp(i,2)
 1007 sd(i,6)=pzbz(i)*sp(i,3)+c(i)

      do 1008 j=1,3
      do 1008 i=1,n
 1008 c1(i,j)=ps0(i,j)/twozeta(i)

      do 1009 i=1,n
      fs0(i,01)=pxax(i)*ds0(i,1)+c1(i,1)+c1(i,1)
      fs0(i,02)=pyay(i)*ds0(i,1)
      fs0(i,03)=pzaz(i)*ds0(i,1)
      fs0(i,04)=pyay(i)*ds0(i,2)+c1(i,1)
      fs0(i,05)=pzaz(i)*ds0(i,2)
      fs0(i,06)=pzaz(i)*ds0(i,3)+c1(i,1)
      fs0(i,07)=pyay(i)*ds0(i,4)+c1(i,2)+c1(i,2)
      fs0(i,08)=pzaz(i)*ds0(i,4)
      fs0(i,09)=pzaz(i)*ds0(i,5)+c1(i,2)
 1009 fs0(i,10)=pzaz(i)*ds0(i,6)+c1(i,3)+c1(i,3)

      do 1010 i=1,n
      dp(i,01)=pxbx(i)*ds0(i,1)+c1(i,1)*two
      dp(i,02)=pyby(i)*ds0(i,1)
      dp(i,03)=pzbz(i)*ds0(i,1)
      dp(i,04)=pxbx(i)*ds0(i,2)+c1(i,2)
      dp(i,05)=pyby(i)*ds0(i,2)+c1(i,1)
      dp(i,06)=pzbz(i)*ds0(i,2)
      dp(i,07)=pxbx(i)*ds0(i,3)+c1(i,3)
      dp(i,08)=pyby(i)*ds0(i,3)
      dp(i,09)=pzbz(i)*ds0(i,3)+c1(i,1)
      dp(i,10)=pxbx(i)*ds0(i,4)
      dp(i,11)=pyby(i)*ds0(i,4)+c1(i,2)*two
      dp(i,12)=pzbz(i)*ds0(i,4)
      dp(i,13)=pxbx(i)*ds0(i,5)
      dp(i,14)=pyby(i)*ds0(i,5)+c1(i,3)
      dp(i,15)=pzbz(i)*ds0(i,5)+c1(i,2)
      dp(i,16)=pxbx(i)*ds0(i,6)
      dp(i,17)=pyby(i)*ds0(i,6)
 1010 dp(i,18)=pzbz(i)*ds0(i,6)+c1(i,3)*two

      do 1011 j=1,3
      do 1011 i=1,n
 1011 c1(i,j)=sp(i,j)/twozeta(i)

      do 1012 i=1,n
      pd(i,01)=pxax(i)*sd(i,1)+c1(i,1)+c1(i,1)
      pd(i,02)=pxax(i)*sd(i,2)+c1(i,2)
      pd(i,03)=pxax(i)*sd(i,3)+c1(i,3)
      pd(i,04)=pxax(i)*sd(i,4)
      pd(i,05)=pxax(i)*sd(i,5)
      pd(i,06)=pxax(i)*sd(i,6)
      pd(i,07)=pyay(i)*sd(i,1)
      pd(i,08)=pyay(i)*sd(i,2)+c1(i,1)
      pd(i,09)=pyay(i)*sd(i,3)
      pd(i,10)=pyay(i)*sd(i,4)+c1(i,2)+c1(i,2)
      pd(i,11)=pyay(i)*sd(i,5)+c1(i,3)
      pd(i,12)=pyay(i)*sd(i,6)
      pd(i,13)=pzaz(i)*sd(i,1)
      pd(i,14)=pzaz(i)*sd(i,2)
      pd(i,15)=pzaz(i)*sd(i,3)+c1(i,1)
      pd(i,16)=pzaz(i)*sd(i,4)
      pd(i,17)=pzaz(i)*sd(i,5)+c1(i,2)
 1012 pd(i,18)=pzaz(i)*sd(i,6)+c1(i,3)+c1(i,3)

      do 1013 j=1,6
      do 1013 i=1,n
 1013 c1(i,j)=ds0(i,j)/twozeta(i)

      do 1014 i=1,n
      fp(i,01)=pxbx(i)*fs0(i,01)+c1(i,1)*three
      fp(i,02)=pyby(i)*fs0(i,01)
      fp(i,03)=pzbz(i)*fs0(i,01)
      fp(i,04)=pxbx(i)*fs0(i,02)+c1(i,2)*two
      fp(i,05)=pyby(i)*fs0(i,02)+c1(i,1)
      fp(i,06)=pzbz(i)*fs0(i,02)
      fp(i,07)=pxbx(i)*fs0(i,03)+c1(i,3)*two
      fp(i,08)=pyby(i)*fs0(i,03)
      fp(i,09)=pzbz(i)*fs0(i,03)+c1(i,1)
      fp(i,10)=pxbx(i)*fs0(i,04)+c1(i,4)
      fp(i,11)=pyby(i)*fs0(i,04)+c1(i,2)*two
      fp(i,12)=pzbz(i)*fs0(i,04)
      fp(i,13)=pxbx(i)*fs0(i,05)+c1(i,5)
      fp(i,14)=pyby(i)*fs0(i,05)+c1(i,3)
      fp(i,15)=pzbz(i)*fs0(i,05)+c1(i,2)
      fp(i,16)=pxbx(i)*fs0(i,06)+c1(i,6)
      fp(i,17)=pyby(i)*fs0(i,06)
      fp(i,18)=pzbz(i)*fs0(i,06)+c1(i,3)*two
      fp(i,19)=pxbx(i)*fs0(i,07)
      fp(i,20)=pyby(i)*fs0(i,07)+c1(i,4)*three
      fp(i,21)=pzbz(i)*fs0(i,07)
      fp(i,22)=pxbx(i)*fs0(i,08)
      fp(i,23)=pyby(i)*fs0(i,08)+c1(i,5)*two
      fp(i,24)=pzbz(i)*fs0(i,08)+c1(i,4)
      fp(i,25)=pxbx(i)*fs0(i,09)
      fp(i,26)=pyby(i)*fs0(i,09)+c1(i,6)
      fp(i,27)=pzbz(i)*fs0(i,09)+c1(i,5)*two
      fp(i,28)=pxbx(i)*fs0(i,10)
      fp(i,29)=pyby(i)*fs0(i,10)
 1014 fp(i,30)=pzbz(i)*fs0(i,10)+c1(i,6)*three

      do 1015 j=1,10
      do 1015 i=1,n
 1015 c1(i,j)=fs0(i,j)/twozeta(i)

      do 1016 j=1,18
      do 1016 i=1,n
 1016 c2(i,j)=dp(i,j)/twozeta(i)

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

      do 1018 j=1,18
      do 1018 i=1,m
 1018 pdcond(i,j)=zero

      do 1019 j=1,60
      do 1019 i=1,m
 1019 fdcond(i,j)=zero

      do 1020 j=1,18
      do 1020 i=1,n
 1020 pdcond(icontr(i),j)=pdcond(icontr(i),j)+pd(i,j)

      do 1021 j=1,60
      do 1021 i=1,n
 1021 fdcond(icontr(i),j)=fdcond(icontr(i),j)+fd(i,j)*twozetaa(i)

      do 1022 i=1,m
      overlap(inddyy(i,01),1)=overlap(inddyy(i,01),1)-pdcond(i,01)*two
      overlap(inddyy(i,07),1)=overlap(inddyy(i,07),1)-pdcond(i,07)
      overlap(inddyy(i,08),1)=overlap(inddyy(i,08),1)-pdcond(i,08)
      overlap(inddyy(i,13),1)=overlap(inddyy(i,13),1)-pdcond(i,13)
      overlap(inddyy(i,14),1)=overlap(inddyy(i,14),1)-pdcond(i,14)
      overlap(inddyy(i,15),1)=overlap(inddyy(i,15),1)-pdcond(i,15)
      overlap(inddyy(i,07),2)=overlap(inddyy(i,07),2)-pdcond(i,01)
      overlap(inddyy(i,08),2)=overlap(inddyy(i,08),2)-pdcond(i,02)
      overlap(inddyy(i,19),2)=overlap(inddyy(i,19),2)-pdcond(i,07)*two
      overlap(inddyy(i,20),2)=overlap(inddyy(i,20),2)-pdcond(i,08)*two
      overlap(inddyy(i,21),2)=overlap(inddyy(i,21),2)-pdcond(i,09)*two
      overlap(inddyy(i,22),2)=overlap(inddyy(i,22),2)-pdcond(i,10)*two
      overlap(inddyy(i,25),2)=overlap(inddyy(i,25),2)-pdcond(i,13)
      overlap(inddyy(i,26),2)=overlap(inddyy(i,26),2)-pdcond(i,14)
      overlap(inddyy(i,27),2)=overlap(inddyy(i,27),2)-pdcond(i,15)
      overlap(inddyy(i,28),2)=overlap(inddyy(i,28),2)-pdcond(i,16)
      overlap(inddyy(i,29),2)=overlap(inddyy(i,29),2)-pdcond(i,17)
      overlap(inddyy(i,13),3)=overlap(inddyy(i,13),3)-pdcond(i,01)
      overlap(inddyy(i,14),3)=overlap(inddyy(i,14),3)-pdcond(i,02)
      overlap(inddyy(i,15),3)=overlap(inddyy(i,15),3)-pdcond(i,03)
      overlap(inddyy(i,25),3)=overlap(inddyy(i,25),3)-pdcond(i,07)
      overlap(inddyy(i,26),3)=overlap(inddyy(i,26),3)-pdcond(i,08)
      overlap(inddyy(i,27),3)=overlap(inddyy(i,27),3)-pdcond(i,09)
      overlap(inddyy(i,28),3)=overlap(inddyy(i,28),3)-pdcond(i,10)
      overlap(inddyy(i,29),3)=overlap(inddyy(i,29),3)-pdcond(i,11)
      overlap(inddyy(i,31),3)=overlap(inddyy(i,31),3)-pdcond(i,13)*two
      overlap(inddyy(i,32),3)=overlap(inddyy(i,32),3)-pdcond(i,14)*two
      overlap(inddyy(i,33),3)=overlap(inddyy(i,33),3)-pdcond(i,15)*two
      overlap(inddyy(i,34),3)=overlap(inddyy(i,34),3)-pdcond(i,16)*two
      overlap(inddyy(i,35),3)=overlap(inddyy(i,35),3)-pdcond(i,17)*two
 1022 overlap(inddyy(i,36),3)=overlap(inddyy(i,36),3)-pdcond(i,18)*two

      do 1023 i=1,m
      if(iminj(i).ne.0) then

        overlap(inddyy(i,02),1)=overlap(inddyy(i,02),1)-pdcond(i,02)*two
        overlap(inddyy(i,03),1)=overlap(inddyy(i,03),1)-pdcond(i,03)*two
        overlap(inddyy(i,04),1)=overlap(inddyy(i,04),1)-pdcond(i,04)*two
        overlap(inddyy(i,05),1)=overlap(inddyy(i,05),1)-pdcond(i,05)*two
        overlap(inddyy(i,06),1)=overlap(inddyy(i,06),1)-pdcond(i,06)*two
        overlap(inddyy(i,09),1)=overlap(inddyy(i,09),1)-pdcond(i,09)
        overlap(inddyy(i,10),1)=overlap(inddyy(i,10),1)-pdcond(i,10)
        overlap(inddyy(i,11),1)=overlap(inddyy(i,11),1)-pdcond(i,11)
        overlap(inddyy(i,12),1)=overlap(inddyy(i,12),1)-pdcond(i,12)
        overlap(inddyy(i,16),1)=overlap(inddyy(i,16),1)-pdcond(i,16)
        overlap(inddyy(i,17),1)=overlap(inddyy(i,17),1)-pdcond(i,17)
        overlap(inddyy(i,18),1)=overlap(inddyy(i,18),1)-pdcond(i,18)
        overlap(inddyy(i,09),2)=overlap(inddyy(i,09),2)-pdcond(i,03)
        overlap(inddyy(i,10),2)=overlap(inddyy(i,10),2)-pdcond(i,04)
        overlap(inddyy(i,11),2)=overlap(inddyy(i,11),2)-pdcond(i,05)
        overlap(inddyy(i,12),2)=overlap(inddyy(i,12),2)-pdcond(i,06)
        overlap(inddyy(i,23),2)=overlap(inddyy(i,23),2)-pdcond(i,11)*two
        overlap(inddyy(i,24),2)=overlap(inddyy(i,24),2)-pdcond(i,12)*two
        overlap(inddyy(i,30),2)=overlap(inddyy(i,30),2)-pdcond(i,18)
        overlap(inddyy(i,16),3)=overlap(inddyy(i,16),3)-pdcond(i,04)
        overlap(inddyy(i,17),3)=overlap(inddyy(i,17),3)-pdcond(i,05)
        overlap(inddyy(i,18),3)=overlap(inddyy(i,18),3)-pdcond(i,06)
        overlap(inddyy(i,30),3)=overlap(inddyy(i,30),3)-pdcond(i,12)

                        endif
 1023 continue

      do 1024 i=1,m
      overlap(inddyy(i,01),1)=overlap(inddyy(i,01),1)+fdcond(i,01)
      overlap(inddyy(i,07),1)=overlap(inddyy(i,07),1)+fdcond(i,07)
      overlap(inddyy(i,08),1)=overlap(inddyy(i,08),1)+fdcond(i,08)
      overlap(inddyy(i,13),1)=overlap(inddyy(i,13),1)+fdcond(i,13)
      overlap(inddyy(i,14),1)=overlap(inddyy(i,14),1)+fdcond(i,14)
      overlap(inddyy(i,15),1)=overlap(inddyy(i,15),1)+fdcond(i,15)
      overlap(inddyy(i,19),1)=overlap(inddyy(i,19),1)+fdcond(i,19)
      overlap(inddyy(i,20),1)=overlap(inddyy(i,20),1)+fdcond(i,20)
      overlap(inddyy(i,21),1)=overlap(inddyy(i,21),1)+fdcond(i,21)
      overlap(inddyy(i,22),1)=overlap(inddyy(i,22),1)+fdcond(i,22)
      overlap(inddyy(i,25),1)=overlap(inddyy(i,25),1)+fdcond(i,25)
      overlap(inddyy(i,26),1)=overlap(inddyy(i,26),1)+fdcond(i,26)
      overlap(inddyy(i,27),1)=overlap(inddyy(i,27),1)+fdcond(i,27)
      overlap(inddyy(i,28),1)=overlap(inddyy(i,28),1)+fdcond(i,28)
      overlap(inddyy(i,29),1)=overlap(inddyy(i,29),1)+fdcond(i,29)
      overlap(inddyy(i,31),1)=overlap(inddyy(i,31),1)+fdcond(i,31)
      overlap(inddyy(i,32),1)=overlap(inddyy(i,32),1)+fdcond(i,32)
      overlap(inddyy(i,33),1)=overlap(inddyy(i,33),1)+fdcond(i,33)
      overlap(inddyy(i,34),1)=overlap(inddyy(i,34),1)+fdcond(i,34)
      overlap(inddyy(i,35),1)=overlap(inddyy(i,35),1)+fdcond(i,35)
      overlap(inddyy(i,36),1)=overlap(inddyy(i,36),1)+fdcond(i,36)
      overlap(inddyy(i,01),2)=overlap(inddyy(i,01),2)+fdcond(i,07)
      overlap(inddyy(i,07),2)=overlap(inddyy(i,07),2)+fdcond(i,19)
      overlap(inddyy(i,08),2)=overlap(inddyy(i,08),2)+fdcond(i,20)
      overlap(inddyy(i,13),2)=overlap(inddyy(i,13),2)+fdcond(i,25)
      overlap(inddyy(i,14),2)=overlap(inddyy(i,14),2)+fdcond(i,26)
      overlap(inddyy(i,15),2)=overlap(inddyy(i,15),2)+fdcond(i,27)
      overlap(inddyy(i,19),2)=overlap(inddyy(i,19),2)+fdcond(i,37)
      overlap(inddyy(i,20),2)=overlap(inddyy(i,20),2)+fdcond(i,38)
      overlap(inddyy(i,21),2)=overlap(inddyy(i,21),2)+fdcond(i,39)
      overlap(inddyy(i,22),2)=overlap(inddyy(i,22),2)+fdcond(i,40)
      overlap(inddyy(i,25),2)=overlap(inddyy(i,25),2)+fdcond(i,43)
      overlap(inddyy(i,26),2)=overlap(inddyy(i,26),2)+fdcond(i,44)
      overlap(inddyy(i,27),2)=overlap(inddyy(i,27),2)+fdcond(i,45)
      overlap(inddyy(i,28),2)=overlap(inddyy(i,28),2)+fdcond(i,46)
      overlap(inddyy(i,29),2)=overlap(inddyy(i,29),2)+fdcond(i,47)
      overlap(inddyy(i,31),2)=overlap(inddyy(i,31),2)+fdcond(i,49)
      overlap(inddyy(i,32),2)=overlap(inddyy(i,32),2)+fdcond(i,50)
      overlap(inddyy(i,33),2)=overlap(inddyy(i,33),2)+fdcond(i,51)
      overlap(inddyy(i,34),2)=overlap(inddyy(i,34),2)+fdcond(i,52)
      overlap(inddyy(i,35),2)=overlap(inddyy(i,35),2)+fdcond(i,53)
      overlap(inddyy(i,36),2)=overlap(inddyy(i,36),2)+fdcond(i,54)
      overlap(inddyy(i,01),3)=overlap(inddyy(i,01),3)+fdcond(i,13)
      overlap(inddyy(i,07),3)=overlap(inddyy(i,07),3)+fdcond(i,25)
      overlap(inddyy(i,08),3)=overlap(inddyy(i,08),3)+fdcond(i,26)
      overlap(inddyy(i,13),3)=overlap(inddyy(i,13),3)+fdcond(i,31)
      overlap(inddyy(i,14),3)=overlap(inddyy(i,14),3)+fdcond(i,32)
      overlap(inddyy(i,15),3)=overlap(inddyy(i,15),3)+fdcond(i,33)
      overlap(inddyy(i,19),3)=overlap(inddyy(i,19),3)+fdcond(i,43)
      overlap(inddyy(i,20),3)=overlap(inddyy(i,20),3)+fdcond(i,44)
      overlap(inddyy(i,21),3)=overlap(inddyy(i,21),3)+fdcond(i,45)
      overlap(inddyy(i,22),3)=overlap(inddyy(i,22),3)+fdcond(i,46)
      overlap(inddyy(i,25),3)=overlap(inddyy(i,25),3)+fdcond(i,49)
      overlap(inddyy(i,26),3)=overlap(inddyy(i,26),3)+fdcond(i,50)
      overlap(inddyy(i,27),3)=overlap(inddyy(i,27),3)+fdcond(i,51)
      overlap(inddyy(i,28),3)=overlap(inddyy(i,28),3)+fdcond(i,52)
      overlap(inddyy(i,29),3)=overlap(inddyy(i,29),3)+fdcond(i,53)
      overlap(inddyy(i,31),3)=overlap(inddyy(i,31),3)+fdcond(i,55)
      overlap(inddyy(i,32),3)=overlap(inddyy(i,32),3)+fdcond(i,56)
      overlap(inddyy(i,33),3)=overlap(inddyy(i,33),3)+fdcond(i,57)
      overlap(inddyy(i,34),3)=overlap(inddyy(i,34),3)+fdcond(i,58)
      overlap(inddyy(i,35),3)=overlap(inddyy(i,35),3)+fdcond(i,59)
 1024 overlap(inddyy(i,36),3)=overlap(inddyy(i,36),3)+fdcond(i,60)


      do 1025 i=1,m
      if(iminj(i).ne.0) then

            overlap(inddyy(i,02),1)=overlap(inddyy(i,02),1)+fdcond(i,02)
            overlap(inddyy(i,03),1)=overlap(inddyy(i,03),1)+fdcond(i,03)
            overlap(inddyy(i,04),1)=overlap(inddyy(i,04),1)+fdcond(i,04)
            overlap(inddyy(i,05),1)=overlap(inddyy(i,05),1)+fdcond(i,05)
            overlap(inddyy(i,06),1)=overlap(inddyy(i,06),1)+fdcond(i,06)
            overlap(inddyy(i,09),1)=overlap(inddyy(i,09),1)+fdcond(i,09)
            overlap(inddyy(i,10),1)=overlap(inddyy(i,10),1)+fdcond(i,10)
            overlap(inddyy(i,11),1)=overlap(inddyy(i,11),1)+fdcond(i,11)
            overlap(inddyy(i,12),1)=overlap(inddyy(i,12),1)+fdcond(i,12)
            overlap(inddyy(i,16),1)=overlap(inddyy(i,16),1)+fdcond(i,16)
            overlap(inddyy(i,17),1)=overlap(inddyy(i,17),1)+fdcond(i,17)
            overlap(inddyy(i,18),1)=overlap(inddyy(i,18),1)+fdcond(i,18)
            overlap(inddyy(i,23),1)=overlap(inddyy(i,23),1)+fdcond(i,23)
            overlap(inddyy(i,24),1)=overlap(inddyy(i,24),1)+fdcond(i,24)
            overlap(inddyy(i,30),1)=overlap(inddyy(i,30),1)+fdcond(i,30)
            overlap(inddyy(i,02),2)=overlap(inddyy(i,02),2)+fdcond(i,08)
            overlap(inddyy(i,03),2)=overlap(inddyy(i,03),2)+fdcond(i,09)
            overlap(inddyy(i,04),2)=overlap(inddyy(i,04),2)+fdcond(i,10)
            overlap(inddyy(i,05),2)=overlap(inddyy(i,05),2)+fdcond(i,11)
            overlap(inddyy(i,06),2)=overlap(inddyy(i,06),2)+fdcond(i,12)
            overlap(inddyy(i,09),2)=overlap(inddyy(i,09),2)+fdcond(i,21)
            overlap(inddyy(i,10),2)=overlap(inddyy(i,10),2)+fdcond(i,22)
            overlap(inddyy(i,11),2)=overlap(inddyy(i,11),2)+fdcond(i,23)
            overlap(inddyy(i,12),2)=overlap(inddyy(i,12),2)+fdcond(i,24)
            overlap(inddyy(i,16),2)=overlap(inddyy(i,16),2)+fdcond(i,28)
            overlap(inddyy(i,17),2)=overlap(inddyy(i,17),2)+fdcond(i,29)
            overlap(inddyy(i,18),2)=overlap(inddyy(i,18),2)+fdcond(i,30)
            overlap(inddyy(i,23),2)=overlap(inddyy(i,23),2)+fdcond(i,41)
            overlap(inddyy(i,24),2)=overlap(inddyy(i,24),2)+fdcond(i,42)
            overlap(inddyy(i,30),2)=overlap(inddyy(i,30),2)+fdcond(i,48)
            overlap(inddyy(i,02),3)=overlap(inddyy(i,02),3)+fdcond(i,14)
            overlap(inddyy(i,03),3)=overlap(inddyy(i,03),3)+fdcond(i,15)
            overlap(inddyy(i,04),3)=overlap(inddyy(i,04),3)+fdcond(i,16)
            overlap(inddyy(i,05),3)=overlap(inddyy(i,05),3)+fdcond(i,17)
            overlap(inddyy(i,06),3)=overlap(inddyy(i,06),3)+fdcond(i,18)
            overlap(inddyy(i,09),3)=overlap(inddyy(i,09),3)+fdcond(i,27)
            overlap(inddyy(i,10),3)=overlap(inddyy(i,10),3)+fdcond(i,28)
            overlap(inddyy(i,11),3)=overlap(inddyy(i,11),3)+fdcond(i,29)
            overlap(inddyy(i,12),3)=overlap(inddyy(i,12),3)+fdcond(i,30)
            overlap(inddyy(i,16),3)=overlap(inddyy(i,16),3)+fdcond(i,34)
            overlap(inddyy(i,17),3)=overlap(inddyy(i,17),3)+fdcond(i,35)
            overlap(inddyy(i,18),3)=overlap(inddyy(i,18),3)+fdcond(i,36)
            overlap(inddyy(i,23),3)=overlap(inddyy(i,23),3)+fdcond(i,47)
            overlap(inddyy(i,24),3)=overlap(inddyy(i,24),3)+fdcond(i,48)
            overlap(inddyy(i,30),3)=overlap(inddyy(i,30),3)+fdcond(i,54)

                        endif
 1025 continue

      do 1026 i=1,n
 1026 ss(i)=eta(i)*(three-two*etarr(i))*ss0(i)

      do 1027 i=1,n
      ps(i,1)=pxax(i)*ss(i)+eta2(i)*ps0(i,1)
      ps(i,2)=pyay(i)*ss(i)+eta2(i)*ps0(i,2)
 1027 ps(i,3)=pzaz(i)*ss(i)+eta2(i)*ps0(i,3)

      do 1028 i=1,n
      sp(i,1)=pxbx(i)*ss(i)+eta2(i)*sp(i,1)
      sp(i,2)=pyby(i)*ss(i)+eta2(i)*sp(i,2)
 1028 sp(i,3)=pzbz(i)*ss(i)+eta2(i)*sp(i,3)

      do 1029 i=1,n
 1029 c(i)=ss(i)/twozeta(i)-eta2(i)*ss0(i)/twozetaa(i)

      do 1030 i=1,n
      ds(i,1)=pxax(i)*ps(i,1)+eta2(i)*ds0(i,1)+c(i)
      ds(i,2)=pyay(i)*ps(i,1)+eta2(i)*ds0(i,2)
      ds(i,3)=pzaz(i)*ps(i,1)+eta2(i)*ds0(i,3)
      ds(i,4)=pyay(i)*ps(i,2)+eta2(i)*ds0(i,4)+c(i)
      ds(i,5)=pzaz(i)*ps(i,2)+eta2(i)*ds0(i,5)
 1030 ds(i,6)=pzaz(i)*ps(i,3)+eta2(i)*ds0(i,6)+c(i)

      do 1031 i=1,n
 1031 c(i)=ss(i)/twozeta(i)-eta2(i)*ss0(i)/twozetab(i)

      do 1032 i=1,n
      sd(i,1)=pxbx(i)*sp(i,1)+eta2(i)*sd(i,1)+c(i)
      sd(i,2)=pyby(i)*sp(i,1)+eta2(i)*sd(i,2)
      sd(i,3)=pzbz(i)*sp(i,1)+eta2(i)*sd(i,3)
      sd(i,4)=pyby(i)*sp(i,2)+eta2(i)*sd(i,4)+c(i)
      sd(i,5)=pzbz(i)*sp(i,2)+eta2(i)*sd(i,5)
 1032 sd(i,6)=pzbz(i)*sp(i,3)+eta2(i)*sd(i,6)+c(i)

      do 1033 j=1,3
      do 1033 i=1,n
 1033 c1(i,j)=ps(i,j)/twozeta(i)-eta2(i)*ps0(i,j)/twozetaa(i)

      do 1034 i=1,n
      fs(i,01)=pxax(i)*ds(i,1)+eta2(i)*fs0(i,01)+c1(i,1)+c1(i,1)
      fs(i,02)=pyay(i)*ds(i,1)+eta2(i)*fs0(i,02)
      fs(i,03)=pzaz(i)*ds(i,1)+eta2(i)*fs0(i,03)
      fs(i,04)=pyay(i)*ds(i,2)+eta2(i)*fs0(i,04)+c1(i,1)
      fs(i,05)=pzaz(i)*ds(i,2)+eta2(i)*fs0(i,05)
      fs(i,06)=pzaz(i)*ds(i,3)+eta2(i)*fs0(i,06)+c1(i,1)
      fs(i,07)=pyay(i)*ds(i,4)+eta2(i)*fs0(i,07)+c1(i,2)+c1(i,2)
      fs(i,08)=pzaz(i)*ds(i,4)+eta2(i)*fs0(i,08)
      fs(i,09)=pzaz(i)*ds(i,5)+eta2(i)*fs0(i,09)+c1(i,2)
 1034 fs(i,10)=pzaz(i)*ds(i,6)+eta2(i)*fs0(i,10)+c1(i,3)+c1(i,3)

      do 1035 j=1,3
      do 1035 i=1,n
 1035 c1(i,j)=ps(i,j)/twozeta(i)

      do 1036 i=1,n
      dp(i,01)=pxbx(i)*ds(i,1)+eta2(i)*dp(i,01)+c1(i,1)*two
      dp(i,02)=pyby(i)*ds(i,1)+eta2(i)*dp(i,02)
      dp(i,03)=pzbz(i)*ds(i,1)+eta2(i)*dp(i,03)
      dp(i,04)=pxbx(i)*ds(i,2)+eta2(i)*dp(i,04)+c1(i,2)
      dp(i,05)=pyby(i)*ds(i,2)+eta2(i)*dp(i,05)+c1(i,1)
      dp(i,06)=pzbz(i)*ds(i,2)+eta2(i)*dp(i,06)
      dp(i,07)=pxbx(i)*ds(i,3)+eta2(i)*dp(i,07)+c1(i,3)
      dp(i,08)=pyby(i)*ds(i,3)+eta2(i)*dp(i,08)
      dp(i,09)=pzbz(i)*ds(i,3)+eta2(i)*dp(i,09)+c1(i,1)
      dp(i,10)=pxbx(i)*ds(i,4)+eta2(i)*dp(i,10)
      dp(i,11)=pyby(i)*ds(i,4)+eta2(i)*dp(i,11)+c1(i,2)*two
      dp(i,12)=pzbz(i)*ds(i,4)+eta2(i)*dp(i,12)
      dp(i,13)=pxbx(i)*ds(i,5)+eta2(i)*dp(i,13)
      dp(i,14)=pyby(i)*ds(i,5)+eta2(i)*dp(i,14)+c1(i,3)
      dp(i,15)=pzbz(i)*ds(i,5)+eta2(i)*dp(i,15)+c1(i,2)
      dp(i,16)=pxbx(i)*ds(i,6)+eta2(i)*dp(i,16)
      dp(i,17)=pyby(i)*ds(i,6)+eta2(i)*dp(i,17)
 1036 dp(i,18)=pzbz(i)*ds(i,6)+eta2(i)*dp(i,18)+c1(i,3)*two

      do 1037 j=1,3
      do 1037 i=1,n
 1037 c1(i,j)=sp(i,j)/twozeta(i)

      do 1038 i=1,n
      pd(i,01)=pxax(i)*sd(i,1)+eta2(i)*pd(i,01)+c1(i,1)+c1(i,1)
      pd(i,02)=pxax(i)*sd(i,2)+eta2(i)*pd(i,02)+c1(i,2)
      pd(i,03)=pxax(i)*sd(i,3)+eta2(i)*pd(i,03)+c1(i,3)
      pd(i,04)=pxax(i)*sd(i,4)+eta2(i)*pd(i,04)
      pd(i,05)=pxax(i)*sd(i,5)+eta2(i)*pd(i,05)
      pd(i,06)=pxax(i)*sd(i,6)+eta2(i)*pd(i,06)
      pd(i,07)=pyay(i)*sd(i,1)+eta2(i)*pd(i,07)
      pd(i,08)=pyay(i)*sd(i,2)+eta2(i)*pd(i,08)+c1(i,1)
      pd(i,09)=pyay(i)*sd(i,3)+eta2(i)*pd(i,09)
      pd(i,10)=pyay(i)*sd(i,4)+eta2(i)*pd(i,10)+c1(i,2)+c1(i,2)
      pd(i,11)=pyay(i)*sd(i,5)+eta2(i)*pd(i,11)+c1(i,3)
      pd(i,12)=pyay(i)*sd(i,6)+eta2(i)*pd(i,12)
      pd(i,13)=pzaz(i)*sd(i,1)+eta2(i)*pd(i,13)
      pd(i,14)=pzaz(i)*sd(i,2)+eta2(i)*pd(i,14)
      pd(i,15)=pzaz(i)*sd(i,3)+eta2(i)*pd(i,15)+c1(i,1)
      pd(i,16)=pzaz(i)*sd(i,4)+eta2(i)*pd(i,16)
      pd(i,17)=pzaz(i)*sd(i,5)+eta2(i)*pd(i,17)+c1(i,2)
 1038 pd(i,18)=pzaz(i)*sd(i,6)+eta2(i)*pd(i,18)+c1(i,3)+c1(i,3)

      do 1039 j=1,6
      do 1039 i=1,n
 1039 c1(i,j)=ds(i,j)/twozeta(i)

      do 1040 i=1,n
      fp(i,01)=pxbx(i)*fs(i,01)+eta2(i)*fp(i,01)+c1(i,1)*three
      fp(i,02)=pyby(i)*fs(i,01)+eta2(i)*fp(i,02)
      fp(i,03)=pzbz(i)*fs(i,01)+eta2(i)*fp(i,03)
      fp(i,04)=pxbx(i)*fs(i,02)+eta2(i)*fp(i,04)+c1(i,2)*two
      fp(i,05)=pyby(i)*fs(i,02)+eta2(i)*fp(i,05)+c1(i,1)
      fp(i,06)=pzbz(i)*fs(i,02)+eta2(i)*fp(i,06)
      fp(i,07)=pxbx(i)*fs(i,03)+eta2(i)*fp(i,07)+c1(i,3)*two
      fp(i,08)=pyby(i)*fs(i,03)+eta2(i)*fp(i,08)
      fp(i,09)=pzbz(i)*fs(i,03)+eta2(i)*fp(i,09)+c1(i,1)
      fp(i,10)=pxbx(i)*fs(i,04)+eta2(i)*fp(i,10)+c1(i,4)
      fp(i,11)=pyby(i)*fs(i,04)+eta2(i)*fp(i,11)+c1(i,2)*two
      fp(i,12)=pzbz(i)*fs(i,04)+eta2(i)*fp(i,12)
      fp(i,13)=pxbx(i)*fs(i,05)+eta2(i)*fp(i,13)+c1(i,5)
      fp(i,14)=pyby(i)*fs(i,05)+eta2(i)*fp(i,14)+c1(i,3)
      fp(i,15)=pzbz(i)*fs(i,05)+eta2(i)*fp(i,15)+c1(i,2)
      fp(i,16)=pxbx(i)*fs(i,06)+eta2(i)*fp(i,16)+c1(i,6)
      fp(i,17)=pyby(i)*fs(i,06)+eta2(i)*fp(i,17)
      fp(i,18)=pzbz(i)*fs(i,06)+eta2(i)*fp(i,18)+c1(i,3)*two
      fp(i,19)=pxbx(i)*fs(i,07)+eta2(i)*fp(i,19)
      fp(i,20)=pyby(i)*fs(i,07)+eta2(i)*fp(i,20)+c1(i,4)*three
      fp(i,21)=pzbz(i)*fs(i,07)+eta2(i)*fp(i,21)
      fp(i,22)=pxbx(i)*fs(i,08)+eta2(i)*fp(i,22)
      fp(i,23)=pyby(i)*fs(i,08)+eta2(i)*fp(i,23)+c1(i,5)*two
      fp(i,24)=pzbz(i)*fs(i,08)+eta2(i)*fp(i,24)+c1(i,4)
      fp(i,25)=pxbx(i)*fs(i,09)+eta2(i)*fp(i,25)
      fp(i,26)=pyby(i)*fs(i,09)+eta2(i)*fp(i,26)+c1(i,6)
      fp(i,27)=pzbz(i)*fs(i,09)+eta2(i)*fp(i,27)+c1(i,5)*two
      fp(i,28)=pxbx(i)*fs(i,10)+eta2(i)*fp(i,28)
      fp(i,29)=pyby(i)*fs(i,10)+eta2(i)*fp(i,29)
 1040 fp(i,30)=pzbz(i)*fs(i,10)+eta2(i)*fp(i,30)+c1(i,6)*three

      do 1041 j=1,10
      do 1041 i=1,n
 1041 c1(i,j)=fs(i,j)/twozeta(i)-eta2(i)*fs0(i,j)/twozetab(i)

      do 1042 j=1,18
      do 1042 i=1,n
 1042 c2(i,j)=dp(i,j)/twozeta(i)

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

      do 1044 j=1,18
      do 1044 i=1,m
 1044 pdcond(i,j)=zero

      do 1045 j=1,60
      do 1045 i=1,m
 1045 fdcond(i,j)=zero

      do 1046 j=1,18
      do 1046 i=1,n
 1046 pdcond(icontr(i),j)=pdcond(icontr(i),j)+pd(i,j)

      do 1047 j=1,60
      do 1047 i=1,n
 1047 fdcond(icontr(i),j)=fdcond(icontr(i),j)+fd(i,j)*twozetaa(i)

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

      do 1049 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
 1049 continue

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

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

      do 1052 i=1,m
      grad(i,001)=fdcond(i,01)*density(inddyy(i,01))
      grad(i,002)=fdcond(i,02)*density(inddyy(i,02))
      grad(i,003)=fdcond(i,03)*density(inddyy(i,03))
      grad(i,004)=fdcond(i,04)*density(inddyy(i,04))
      grad(i,005)=fdcond(i,05)*density(inddyy(i,05))
      grad(i,006)=fdcond(i,06)*density(inddyy(i,06))
      grad(i,007)=fdcond(i,07)*density(inddyy(i,07))
      grad(i,008)=fdcond(i,08)*density(inddyy(i,08))
      grad(i,009)=fdcond(i,09)*density(inddyy(i,09))
      grad(i,010)=fdcond(i,10)*density(inddyy(i,10))
      grad(i,011)=fdcond(i,11)*density(inddyy(i,11))
      grad(i,012)=fdcond(i,12)*density(inddyy(i,12))
      grad(i,013)=fdcond(i,13)*density(inddyy(i,13))
      grad(i,014)=fdcond(i,14)*density(inddyy(i,14))
      grad(i,015)=fdcond(i,15)*density(inddyy(i,15))
      grad(i,016)=fdcond(i,16)*density(inddyy(i,16))
      grad(i,017)=fdcond(i,17)*density(inddyy(i,17))
      grad(i,018)=fdcond(i,18)*density(inddyy(i,18))
      grad(i,019)=fdcond(i,19)*density(inddyy(i,19))
      grad(i,020)=fdcond(i,20)*density(inddyy(i,20))
      grad(i,021)=fdcond(i,21)*density(inddyy(i,21))
      grad(i,022)=fdcond(i,22)*density(inddyy(i,22))
      grad(i,023)=fdcond(i,23)*density(inddyy(i,23))
      grad(i,024)=fdcond(i,24)*density(inddyy(i,24))
      grad(i,025)=fdcond(i,25)*density(inddyy(i,25))
      grad(i,026)=fdcond(i,26)*density(inddyy(i,26))
      grad(i,027)=fdcond(i,27)*density(inddyy(i,27))
      grad(i,028)=fdcond(i,28)*density(inddyy(i,28))
      grad(i,029)=fdcond(i,29)*density(inddyy(i,29))
      grad(i,030)=fdcond(i,30)*density(inddyy(i,30))
      grad(i,031)=fdcond(i,31)*density(inddyy(i,31))
      grad(i,032)=fdcond(i,32)*density(inddyy(i,32))
      grad(i,033)=fdcond(i,33)*density(inddyy(i,33))
      grad(i,034)=fdcond(i,34)*density(inddyy(i,34))
      grad(i,035)=fdcond(i,35)*density(inddyy(i,35))
      grad(i,036)=fdcond(i,36)*density(inddyy(i,36))
      grad(i,037)=fdcond(i,07)*density(inddyy(i,01))
      grad(i,038)=fdcond(i,08)*density(inddyy(i,02))
      grad(i,039)=fdcond(i,09)*density(inddyy(i,03))
      grad(i,040)=fdcond(i,10)*density(inddyy(i,04))
      grad(i,041)=fdcond(i,11)*density(inddyy(i,05))
      grad(i,042)=fdcond(i,12)*density(inddyy(i,06))
      grad(i,043)=fdcond(i,19)*density(inddyy(i,07))
      grad(i,044)=fdcond(i,20)*density(inddyy(i,08))
      grad(i,045)=fdcond(i,21)*density(inddyy(i,09))
      grad(i,046)=fdcond(i,22)*density(inddyy(i,10))
      grad(i,047)=fdcond(i,23)*density(inddyy(i,11))
      grad(i,048)=fdcond(i,24)*density(inddyy(i,12))
      grad(i,049)=fdcond(i,25)*density(inddyy(i,13))
      grad(i,050)=fdcond(i,26)*density(inddyy(i,14))
      grad(i,051)=fdcond(i,27)*density(inddyy(i,15))
      grad(i,052)=fdcond(i,28)*density(inddyy(i,16))
      grad(i,053)=fdcond(i,29)*density(inddyy(i,17))
      grad(i,054)=fdcond(i,30)*density(inddyy(i,18))
      grad(i,055)=fdcond(i,37)*density(inddyy(i,19))
      grad(i,056)=fdcond(i,38)*density(inddyy(i,20))
      grad(i,057)=fdcond(i,39)*density(inddyy(i,21))
      grad(i,058)=fdcond(i,40)*density(inddyy(i,22))
      grad(i,059)=fdcond(i,41)*density(inddyy(i,23))
      grad(i,060)=fdcond(i,42)*density(inddyy(i,24))
      grad(i,061)=fdcond(i,43)*density(inddyy(i,25))
      grad(i,062)=fdcond(i,44)*density(inddyy(i,26))
      grad(i,063)=fdcond(i,45)*density(inddyy(i,27))
      grad(i,064)=fdcond(i,46)*density(inddyy(i,28))
      grad(i,065)=fdcond(i,47)*density(inddyy(i,29))
      grad(i,066)=fdcond(i,48)*density(inddyy(i,30))
      grad(i,067)=fdcond(i,49)*density(inddyy(i,31))
      grad(i,068)=fdcond(i,50)*density(inddyy(i,32))
      grad(i,069)=fdcond(i,51)*density(inddyy(i,33))
      grad(i,070)=fdcond(i,52)*density(inddyy(i,34))
      grad(i,071)=fdcond(i,53)*density(inddyy(i,35))
      grad(i,072)=fdcond(i,54)*density(inddyy(i,36))
      grad(i,073)=fdcond(i,13)*density(inddyy(i,01))
      grad(i,074)=fdcond(i,14)*density(inddyy(i,02))
      grad(i,075)=fdcond(i,15)*density(inddyy(i,03))
      grad(i,076)=fdcond(i,16)*density(inddyy(i,04))
      grad(i,077)=fdcond(i,17)*density(inddyy(i,05))
      grad(i,078)=fdcond(i,18)*density(inddyy(i,06))
      grad(i,079)=fdcond(i,25)*density(inddyy(i,07))
      grad(i,080)=fdcond(i,26)*density(inddyy(i,08))
      grad(i,081)=fdcond(i,27)*density(inddyy(i,09))
      grad(i,082)=fdcond(i,28)*density(inddyy(i,10))
      grad(i,083)=fdcond(i,29)*density(inddyy(i,11))
      grad(i,084)=fdcond(i,30)*density(inddyy(i,12))
      grad(i,085)=fdcond(i,31)*density(inddyy(i,13))
      grad(i,086)=fdcond(i,32)*density(inddyy(i,14))
      grad(i,087)=fdcond(i,33)*density(inddyy(i,15))
      grad(i,088)=fdcond(i,34)*density(inddyy(i,16))
      grad(i,089)=fdcond(i,35)*density(inddyy(i,17))
      grad(i,090)=fdcond(i,36)*density(inddyy(i,18))
      grad(i,091)=fdcond(i,43)*density(inddyy(i,19))
      grad(i,092)=fdcond(i,44)*density(inddyy(i,20))
      grad(i,093)=fdcond(i,45)*density(inddyy(i,21))
      grad(i,094)=fdcond(i,46)*density(inddyy(i,22))
      grad(i,095)=fdcond(i,47)*density(inddyy(i,23))
      grad(i,096)=fdcond(i,48)*density(inddyy(i,24))
      grad(i,097)=fdcond(i,49)*density(inddyy(i,25))
      grad(i,098)=fdcond(i,50)*density(inddyy(i,26))
      grad(i,099)=fdcond(i,51)*density(inddyy(i,27))
      grad(i,100)=fdcond(i,52)*density(inddyy(i,28))
      grad(i,101)=fdcond(i,53)*density(inddyy(i,29))
      grad(i,102)=fdcond(i,54)*density(inddyy(i,30))
      grad(i,103)=fdcond(i,55)*density(inddyy(i,31))
      grad(i,104)=fdcond(i,56)*density(inddyy(i,32))
      grad(i,105)=fdcond(i,57)*density(inddyy(i,33))
      grad(i,106)=fdcond(i,58)*density(inddyy(i,34))
      grad(i,107)=fdcond(i,59)*density(inddyy(i,35))
 1052 grad(i,108)=fdcond(i,60)*density(inddyy(i,36))

      do 1053 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,002)=zero
                              grad(i,003)=zero
                              grad(i,004)=zero
                              grad(i,005)=zero
                              grad(i,006)=zero
                              grad(i,009)=zero
                              grad(i,010)=zero
                              grad(i,011)=zero
                              grad(i,012)=zero
                              grad(i,016)=zero
                              grad(i,017)=zero
                              grad(i,018)=zero
                              grad(i,023)=zero
                              grad(i,024)=zero
                              grad(i,030)=zero
                              grad(i,038)=zero
                              grad(i,039)=zero
                              grad(i,040)=zero
                              grad(i,041)=zero
                              grad(i,042)=zero
                              grad(i,045)=zero
                              grad(i,046)=zero
                              grad(i,047)=zero
                              grad(i,048)=zero
                              grad(i,052)=zero
                              grad(i,053)=zero
                              grad(i,054)=zero
                              grad(i,059)=zero
                              grad(i,060)=zero
                              grad(i,066)=zero
                              grad(i,074)=zero
                              grad(i,075)=zero
                              grad(i,076)=zero
                              grad(i,077)=zero
                              grad(i,078)=zero
                              grad(i,081)=zero
                              grad(i,082)=zero
                              grad(i,083)=zero
                              grad(i,084)=zero
                              grad(i,088)=zero
                              grad(i,089)=zero
                              grad(i,090)=zero
                              grad(i,095)=zero
                              grad(i,096)=zero
                              grad(i,102)=zero
                        endif
 1053 continue

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

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

      return
      end
