      subroutine gradovss(ndim,m,n,icontr,inddyy,icentr,density,overlap,
     &                    gradient,zeta,twozetaa,eta,eta2,base,etarr,
     &                    pxax,pyay,pzaz,ss,ps,pscond,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 density(*),overlap(ndim,*),gradient(3,*),zeta(*),
     &          twozetaa(*),eta(*),eta2(*),base(*),etarr(*),pxax(*),
     &          pyay(*),pzaz(*),ss(*),ps(n,*),pscond(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 ss(i)=((pi/zeta(i))**thrhlf)*base(i)

      do 1002 i=1,n
      ps(i,1)=pxax(i)*ss(i)
      ps(i,2)=pyay(i)*ss(i)
 1002 ps(i,3)=pzaz(i)*ss(i)

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

      do 1004 j=1,3
      do 1004 i=1,n
 1004 pscond(icontr(i),j)=pscond(icontr(i),j)+ps(i,j)*twozetaa(i)

      do 1005 j=1,3
      do 1005 i=1,m
 1005 overlap(inddyy(i),j)=overlap(inddyy(i),j)+pscond(i,j)

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

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

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

      do 1009 j=1,3
      do 1009 i=1,n
 1009 pscond(icontr(i),j)=pscond(icontr(i),j)+ps(i,j)*twozetaa(i)

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

      do 1011 j=1,3
      do 1011 i=1,m
 1011 gradient(j,icentr(i,1))=gradient(j,icentr(i,1))+grad(i,j)

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

      return
      end
