      subroutine gradnadp(m,n,ncentr,icontr,inddyy,icentr,g,coord,
     &                    charge,density,gradient,zeta,twozetaa,
     &                    twozetab,px,py,pz,base,pxax,pyay,pzaz,axbx,
     &                    ayby,azbz,pxcx,pycy,pzcz,twozeta,u,ss0,ss1,
     &                    ss2,ss3,ss4,ps0,ps1,ps2,ps3,ds0,ds1,ds2,fs0,
     &                    fs1,gs0,c,c1,pscond,dscond,ppcond,fscond,
     &                    dpcond,gscond,fpcond,ddcond,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 icontr(*),inddyy(m,*),icentr(m,*)

      dimension g(*),coord(3,*),charge(*),density(*),gradient(3,*),
     &          zeta(*),twozetaa(*),twozetab(*),px(*),py(*),pz(*),
     &          base(*),pxax(*),pyay(*),pzaz(*),axbx(*),ayby(*),azbz(*),
     &          pxcx(*),pycy(*),pzcz(*),twozeta(*),u(*),ss0(*),ss1(*),
     &          ss2(*),ss3(*),ss4(*),ps0(n,*),ps1(n,*),ps2(n,*),
     &          ps3(n,*),ds0(n,*),ds1(n,*),ds2(n,*),fs0(n,*),fs1(n,*),
     &          gs0(n,*),c(*),c1(n,*),pscond(m,*),dscond(m,*),
     &          ppcond(m,*),fscond(m,*),dpcond(m,*),gscond(m,*),
     &          fpcond(m,*),ddcond(m,*),grad(m,*)

      data zero,two,three,four,pi/0.0d0,2.0d0,
     &                            3.0d0,4.0d0,3.1415926535898d0/

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

      do 1002 nn=1,ncentr

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

      twopichrg=two*pi*charge(nn)

      do 1003 i=1,n
      pxcx(i)=px(i)-ccx
      pycy(i)=py(i)-ccy
      pzcz(i)=pz(i)-ccz
 1003 u(i)=zeta(i)*(pxcx(i)*pxcx(i)+pycy(i)*pycy(i)+pzcz(i)*pzcz(i))

      do 1004 i=1,n
      call augg(5,u(i),g)
      ss0(i)=base(i)*g(1)*twopichrg/zeta(i)
      ss1(i)=base(i)*g(2)*twopichrg/zeta(i)
      ss2(i)=base(i)*g(3)*twopichrg/zeta(i)
      ss3(i)=base(i)*g(4)*twopichrg/zeta(i)
 1004 ss4(i)=base(i)*g(5)*twopichrg/zeta(i)

      do 1005 i=1,n
      ps0(i,1)=pxax(i)*ss0(i)-pxcx(i)*ss1(i)
      ps0(i,2)=pyay(i)*ss0(i)-pycy(i)*ss1(i)
 1005 ps0(i,3)=pzaz(i)*ss0(i)-pzcz(i)*ss1(i)

      do 1006 i=1,n
      ps1(i,1)=pxax(i)*ss1(i)-pxcx(i)*ss2(i)
      ps1(i,2)=pyay(i)*ss1(i)-pycy(i)*ss2(i)
 1006 ps1(i,3)=pzaz(i)*ss1(i)-pzcz(i)*ss2(i)

      do 1007 i=1,n
      ps2(i,1)=pxax(i)*ss2(i)-pxcx(i)*ss3(i)
      ps2(i,2)=pyay(i)*ss2(i)-pycy(i)*ss3(i)
 1007 ps2(i,3)=pzaz(i)*ss2(i)-pzcz(i)*ss3(i)

      do 1008 i=1,n
      ps3(i,1)=pxax(i)*ss3(i)-pxcx(i)*ss4(i)
      ps3(i,2)=pyay(i)*ss3(i)-pycy(i)*ss4(i)
 1008 ps3(i,3)=pzaz(i)*ss3(i)-pzcz(i)*ss4(i)

      do 1009 i=1,n
 1009 c(i)=(ss0(i)-ss1(i))/twozeta(i)

      do 1010 i=1,n
      ds0(i,1)=pxax(i)*ps0(i,1)-pxcx(i)*ps1(i,1)+c(i)
      ds0(i,2)=pyay(i)*ps0(i,1)-pycy(i)*ps1(i,1)
      ds0(i,3)=pzaz(i)*ps0(i,1)-pzcz(i)*ps1(i,1)
      ds0(i,4)=pyay(i)*ps0(i,2)-pycy(i)*ps1(i,2)+c(i)
      ds0(i,5)=pzaz(i)*ps0(i,2)-pzcz(i)*ps1(i,2)
 1010 ds0(i,6)=pzaz(i)*ps0(i,3)-pzcz(i)*ps1(i,3)+c(i)

      do 1011 i=1,n
 1011 c(i)=(ss1(i)-ss2(i))/twozeta(i)

      do 1012 i=1,n
      ds1(i,1)=pxax(i)*ps1(i,1)-pxcx(i)*ps2(i,1)+c(i)
      ds1(i,2)=pyay(i)*ps1(i,1)-pycy(i)*ps2(i,1)
      ds1(i,3)=pzaz(i)*ps1(i,1)-pzcz(i)*ps2(i,1)
      ds1(i,4)=pyay(i)*ps1(i,2)-pycy(i)*ps2(i,2)+c(i)
      ds1(i,5)=pzaz(i)*ps1(i,2)-pzcz(i)*ps2(i,2)
 1012 ds1(i,6)=pzaz(i)*ps1(i,3)-pzcz(i)*ps2(i,3)+c(i)

      do 1013 i=1,n
 1013 c(i)=(ss2(i)-ss3(i))/twozeta(i)

      do 1014 i=1,n
      ds2(i,1)=pxax(i)*ps2(i,1)-pxcx(i)*ps3(i,1)+c(i)
      ds2(i,2)=pyay(i)*ps2(i,1)-pycy(i)*ps3(i,1)
      ds2(i,3)=pzaz(i)*ps2(i,1)-pzcz(i)*ps3(i,1)
      ds2(i,4)=pyay(i)*ps2(i,2)-pycy(i)*ps3(i,2)+c(i)
      ds2(i,5)=pzaz(i)*ps2(i,2)-pzcz(i)*ps3(i,2)
 1014 ds2(i,6)=pzaz(i)*ps2(i,3)-pzcz(i)*ps3(i,3)+c(i)

      do 1015 j=1,3
      do 1015 i=1,n
 1015 c1(i,j)=(ps0(i,j)-ps1(i,j))/twozeta(i)

      do 1016 i=1,n
      fs0(i,01)=pxax(i)*ds0(i,1)-pxcx(i)*ds1(i,1)+c1(i,1)*two
      fs0(i,02)=pyay(i)*ds0(i,1)-pycy(i)*ds1(i,1)
      fs0(i,03)=pzaz(i)*ds0(i,1)-pzcz(i)*ds1(i,1)
      fs0(i,04)=pyay(i)*ds0(i,2)-pycy(i)*ds1(i,2)+c1(i,1)
      fs0(i,05)=pzaz(i)*ds0(i,2)-pzcz(i)*ds1(i,2)
      fs0(i,06)=pzaz(i)*ds0(i,3)-pzcz(i)*ds1(i,3)+c1(i,1)
      fs0(i,07)=pyay(i)*ds0(i,4)-pycy(i)*ds1(i,4)+c1(i,2)*two
      fs0(i,08)=pzaz(i)*ds0(i,4)-pzcz(i)*ds1(i,4)
      fs0(i,09)=pzaz(i)*ds0(i,5)-pzcz(i)*ds1(i,5)+c1(i,2)
 1016 fs0(i,10)=pzaz(i)*ds0(i,6)-pzcz(i)*ds1(i,6)+c1(i,3)*two

      do 1017 j=1,3
      do 1017 i=1,n
 1017 c1(i,j)=(ps1(i,j)-ps2(i,j))/twozeta(i)

      do 1018 i=1,n
      fs1(i,01)=pxax(i)*ds1(i,1)-pxcx(i)*ds2(i,1)+c1(i,1)*two
      fs1(i,02)=pyay(i)*ds1(i,1)-pycy(i)*ds2(i,1)
      fs1(i,03)=pzaz(i)*ds1(i,1)-pzcz(i)*ds2(i,1)
      fs1(i,04)=pyay(i)*ds1(i,2)-pycy(i)*ds2(i,2)+c1(i,1)
      fs1(i,05)=pzaz(i)*ds1(i,2)-pzcz(i)*ds2(i,2)
      fs1(i,06)=pzaz(i)*ds1(i,3)-pzcz(i)*ds2(i,3)+c1(i,1)
      fs1(i,07)=pyay(i)*ds1(i,4)-pycy(i)*ds2(i,4)+c1(i,2)*two
      fs1(i,08)=pzaz(i)*ds1(i,4)-pzcz(i)*ds2(i,4)
      fs1(i,09)=pzaz(i)*ds1(i,5)-pzcz(i)*ds2(i,5)+c1(i,2)
 1018 fs1(i,10)=pzaz(i)*ds1(i,6)-pzcz(i)*ds2(i,6)+c1(i,3)*two

      do 1019 j=1,6
      do 1019 i=1,n
 1019 c1(i,j)=(ds0(i,j)-ds1(i,j))/twozeta(i)

      do 1020 i=1,n
      gs0(i,01)=pxax(i)*fs0(i,01)-pxcx(i)*fs1(i,01)+c1(i,1)*three
      gs0(i,02)=pyay(i)*fs0(i,01)-pycy(i)*fs1(i,01)
      gs0(i,03)=pzaz(i)*fs0(i,01)-pzcz(i)*fs1(i,01)
      gs0(i,04)=pyay(i)*fs0(i,02)-pycy(i)*fs1(i,02)+c1(i,1)
      gs0(i,05)=pzaz(i)*fs0(i,02)-pzcz(i)*fs1(i,02)
      gs0(i,06)=pzaz(i)*fs0(i,03)-pzcz(i)*fs1(i,03)+c1(i,1)
      gs0(i,07)=pyay(i)*fs0(i,04)-pycy(i)*fs1(i,04)+c1(i,2)*two
      gs0(i,08)=pzaz(i)*fs0(i,04)-pzcz(i)*fs1(i,04)
      gs0(i,09)=pzaz(i)*fs0(i,05)-pzcz(i)*fs1(i,05)+c1(i,2)
      gs0(i,10)=pzaz(i)*fs0(i,06)-pzcz(i)*fs1(i,06)+c1(i,3)*two
      gs0(i,11)=pyay(i)*fs0(i,07)-pycy(i)*fs1(i,07)+c1(i,4)*three
      gs0(i,12)=pzaz(i)*fs0(i,07)-pzcz(i)*fs1(i,07)
      gs0(i,13)=pzaz(i)*fs0(i,08)-pzcz(i)*fs1(i,08)+c1(i,4)
      gs0(i,14)=pzaz(i)*fs0(i,09)-pzcz(i)*fs1(i,09)+c1(i,5)*two
 1020 gs0(i,15)=pzaz(i)*fs0(i,10)-pzcz(i)*fs1(i,10)+c1(i,6)*three

      do 1021 j=1,3
      do 1021 i=1,m
 1021 pscond(i,j)=zero

      do 1022 j=1,6
      do 1022 i=1,m
 1022 dscond(i,j)=zero

      do 1023 j=1,3
      do 1023 i=1,n
 1023 pscond(icontr(i),j)=pscond(icontr(i),j)-ps0(i,j)

      do 1024 j=1,6
      do 1024 i=1,n
 1024 dscond(icontr(i),j)=dscond(icontr(i),j)-ds0(i,j)

      do 1025 i=1,m
      ppcond(i,1)=axbx(i)*pscond(i,1)+dscond(i,1)
      ppcond(i,2)=ayby(i)*pscond(i,1)+dscond(i,2)
      ppcond(i,3)=azbz(i)*pscond(i,1)+dscond(i,3)
      ppcond(i,4)=axbx(i)*pscond(i,2)+dscond(i,2)
      ppcond(i,5)=ayby(i)*pscond(i,2)+dscond(i,4)
      ppcond(i,6)=azbz(i)*pscond(i,2)+dscond(i,5)
      ppcond(i,7)=axbx(i)*pscond(i,3)+dscond(i,3)
      ppcond(i,8)=ayby(i)*pscond(i,3)+dscond(i,5)
 1025 ppcond(i,9)=azbz(i)*pscond(i,3)+dscond(i,6)

      do 1026 i=1,m
      grad(i,01)=ppcond(i,1)*density(inddyy(i,01))*two
      grad(i,02)=ppcond(i,2)*density(inddyy(i,02))*two
      grad(i,03)=ppcond(i,3)*density(inddyy(i,03))*two
      grad(i,04)=ppcond(i,4)*density(inddyy(i,04))
      grad(i,05)=ppcond(i,5)*density(inddyy(i,05))
      grad(i,06)=ppcond(i,6)*density(inddyy(i,06))
      grad(i,07)=ppcond(i,7)*density(inddyy(i,07))
      grad(i,08)=ppcond(i,8)*density(inddyy(i,08))
      grad(i,09)=ppcond(i,9)*density(inddyy(i,09))
      grad(i,10)=ppcond(i,1)*density(inddyy(i,04))
      grad(i,11)=ppcond(i,2)*density(inddyy(i,05))
      grad(i,12)=ppcond(i,3)*density(inddyy(i,06))
      grad(i,13)=ppcond(i,4)*density(inddyy(i,10))*two
      grad(i,14)=ppcond(i,5)*density(inddyy(i,11))*two
      grad(i,15)=ppcond(i,6)*density(inddyy(i,12))*two
      grad(i,16)=ppcond(i,7)*density(inddyy(i,13))
      grad(i,17)=ppcond(i,8)*density(inddyy(i,14))
      grad(i,18)=ppcond(i,9)*density(inddyy(i,15))
      grad(i,19)=ppcond(i,1)*density(inddyy(i,07))
      grad(i,20)=ppcond(i,2)*density(inddyy(i,08))
      grad(i,21)=ppcond(i,3)*density(inddyy(i,09))
      grad(i,22)=ppcond(i,4)*density(inddyy(i,13))
      grad(i,23)=ppcond(i,5)*density(inddyy(i,14))
      grad(i,24)=ppcond(i,6)*density(inddyy(i,15))
      grad(i,25)=ppcond(i,7)*density(inddyy(i,16))*two
      grad(i,26)=ppcond(i,8)*density(inddyy(i,17))*two
 1026 grad(i,27)=ppcond(i,9)*density(inddyy(i,18))*two

      do 1027 k=1,3
      do 1027 j=1,9
      do 1027 i=1,m
 1027 gradient(k,nn)=gradient(k,nn)+grad(i,9*(k-1)+j)

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

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

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

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

      do 1032 j=1,10
      do 1032 i=1,m
 1032 fscond(i,j)=zero

      do 1033 j=1,15
      do 1033 i=1,m
 1033 gscond(i,j)=zero

      do 1034 j=1,10
      do 1034 i=1,n
 1034 fscond(icontr(i),j)=fscond(icontr(i),j)+fs0(i,j)*twozetaa(i)

      do 1035 j=1,15
      do 1035 i=1,n
 1035 gscond(icontr(i),j)=gscond(icontr(i),j)+gs0(i,j)*twozetaa(i)

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

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

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

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

      do 1040 j=1,6
      do 1040 i=1,m
 1040 dscond(i,j)=zero

      do 1041 j=1,10
      do 1041 i=1,m
 1041 fscond(i,j)=zero

      do 1042 j=1,15
      do 1042 i=1,m
 1042 gscond(i,j)=zero

      do 1043 j=1,6
      do 1043 i=1,n
 1043 dscond(icontr(i),j)=dscond(icontr(i),j)+ds0(i,j)*twozetab(i)

      do 1044 j=1,10
      do 1044 i=1,n
 1044 fscond(icontr(i),j)=fscond(icontr(i),j)+fs0(i,j)*twozetab(i)

      do 1045 j=1,15
      do 1045 i=1,n
 1045 gscond(icontr(i),j)=gscond(icontr(i),j)+gs0(i,j)*twozetab(i)

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

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

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

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

      do 1050 k=1,3
      do 1050 j=1,18
      do 1050 i=1,m
 1050 gradient(k,nn)=gradient(k,nn)+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)

 1002 continue

      return
      end
