      subroutine gradnaps(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,ps0,ps1,ds0,c,sscond,pscond,dscond,ppcond,
     &                    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(*),ps0(n,*),ps1(n,*),ds0(n,*),c(*),sscond(*),
     &          pscond(m,*),dscond(m,*),ppcond(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(3,u(i),g)
      ss0(i)=base(i)*g(1)*twopichrg/zeta(i)
      ss1(i)=base(i)*g(2)*twopichrg/zeta(i)
 1004 ss2(i)=base(i)*g(3)*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
 1007 c(i)=(ss0(i)-ss1(i))/twozeta(i)

      do 1008 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)
 1008 ds0(i,6)=pzaz(i)*ps0(i,3)-pzcz(i)*ps1(i,3)+c(i)

      do 1009 i=1,m
 1009 sscond(i)=zero

      do 1010 i=1,n
 1010 sscond(icontr(i))=sscond(icontr(i))-ss0(i)

      do 1011 j=1,3
      do 1011 i=1,m
 1011 grad(i,j)=sscond(i)*density(inddyy(i,j))

      do 1012 j=1,3
      do 1012 i=1,m
 1012 gradient(j,nn)=gradient(j,nn)+grad(i,j)

      do 1013 j=1,3
      do 1013 i=1,m
 1013 gradient(j,icentr(i,1))=gradient(j,icentr(i,1))-grad(i,j)

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

      do 1015 j=1,6
      do 1015 i=1,n
 1015 dscond(icontr(i),j)=dscond(icontr(i),j)+ds0(i,j)*twozetaa(i)

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

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

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

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

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

      do 1021 j=1,3
      do 1021 i=1,n
 1021 pscond(icontr(i),j)=pscond(icontr(i),j)+ps0(i,j)*twozetab(i)

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

      do 1023 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)
 1023 ppcond(i,9)=azbz(i)*pscond(i,3)+dscond(i,6)

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

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

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

 1002 continue

      return
      end
