      subroutine gradnass(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,u,ss0,ss1,ps0,sscond,
     &                    pscond,spcond,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(*),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(*),u(*),ss0(*),ss1(*),ps0(n,*),
     &          sscond(*),pscond(m,*),spcond(m,*),grad(m,*)

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

      do 1001 nn=1,ncentr

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

      twopichrg=two*pi*charge(nn)

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

      do 1003 i=1,n
      call augg(2,u(i),g)
      ss0(i)=base(i)*g(1)*twopichrg/zeta(i)
 1003 ss1(i)=base(i)*g(2)*twopichrg/zeta(i)

      do 1004 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)
 1004 ps0(i,3)=pzaz(i)*ss0(i)-pzcz(i)*ss1(i)

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

      do 1006 j=1,3
      do 1006 i=1,n
 1006 pscond(icontr(i),j)=pscond(icontr(i),j)+ps0(i,j)*twozetaa(i)

      do 1007 j=1,3
      do 1007 i=1,m
 1007 grad(i,j)=pscond(i,j)*density(inddyy(i))

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

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

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

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

      do 1012 i=1,n
 1012 sscond(icontr(i))=sscond(icontr(i))+ss0(i)*twozetab(i)

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

      do 1014 i=1,m
      spcond(i,1)=axbx(i)*sscond(i)+pscond(i,1)
      spcond(i,2)=ayby(i)*sscond(i)+pscond(i,2)
 1014 spcond(i,3)=azbz(i)*sscond(i)+pscond(i,3)

      do 1015 j=1,3
      do 1015 i=1,m
 1015 grad(i,j)=spcond(i,j)*density(inddyy(i))

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

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

 1001 continue

      return
      end
