      subroutine gradovds(ndim,m,n,icontr,inddyy,icentr,density,overlap,
     &                    gradient,zeta,twozetaa,eta,eta2,base,etarr,
     &                    pxax,pyay,pzaz,twozeta,ss0,ss,ps0,ps,ds,fs,c,
     &                    c1,pscond,fscond,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 density(*),overlap(ndim,*),gradient(3,*),zeta(*),
     &          twozetaa(*),eta(*),eta2(*),base(*),etarr(*),pxax(*),
     &          pyay(*),pzaz(*),twozeta(*),ss0(*),ss(*),ps0(n,*),
     &          ps(n,*),ds(n,*),fs(n,*),c(*),c1(n,*),pscond(m,*),
     &          fscond(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
 1004 c(i)=ss0(i)/twozeta(i)

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

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

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

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

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

      do 1010 j=1,3
      do 1010 i=1,n
 1010 pscond(icontr(i),j)=pscond(icontr(i),j)+ps0(i,j)

      do 1011 j=1,10
      do 1011 i=1,n
 1011 fscond(icontr(i),j)=fscond(icontr(i),j)+fs(i,j)*twozetaa(i)

      do 1012 i=1,m
      overlap(inddyy(i,1),1)=overlap(inddyy(i,1),1)-pscond(i,1)*two
      overlap(inddyy(i,2),1)=overlap(inddyy(i,2),1)-pscond(i,2)
      overlap(inddyy(i,3),1)=overlap(inddyy(i,3),1)-pscond(i,3)
      overlap(inddyy(i,2),2)=overlap(inddyy(i,2),2)-pscond(i,1)
      overlap(inddyy(i,4),2)=overlap(inddyy(i,4),2)-pscond(i,2)*two
      overlap(inddyy(i,5),2)=overlap(inddyy(i,5),2)-pscond(i,3)
      overlap(inddyy(i,3),3)=overlap(inddyy(i,3),3)-pscond(i,1)
      overlap(inddyy(i,5),3)=overlap(inddyy(i,5),3)-pscond(i,2)
 1012 overlap(inddyy(i,6),3)=overlap(inddyy(i,6),3)-pscond(i,3)*two

      do 1013 i=1,m
      overlap(inddyy(i,1),1)=overlap(inddyy(i,1),1)+fscond(i,01)
      overlap(inddyy(i,2),1)=overlap(inddyy(i,2),1)+fscond(i,02)
      overlap(inddyy(i,3),1)=overlap(inddyy(i,3),1)+fscond(i,03)
      overlap(inddyy(i,4),1)=overlap(inddyy(i,4),1)+fscond(i,04)
      overlap(inddyy(i,5),1)=overlap(inddyy(i,5),1)+fscond(i,05)
      overlap(inddyy(i,6),1)=overlap(inddyy(i,6),1)+fscond(i,06)
      overlap(inddyy(i,1),2)=overlap(inddyy(i,1),2)+fscond(i,02)
      overlap(inddyy(i,2),2)=overlap(inddyy(i,2),2)+fscond(i,04)
      overlap(inddyy(i,3),2)=overlap(inddyy(i,3),2)+fscond(i,05)
      overlap(inddyy(i,4),2)=overlap(inddyy(i,4),2)+fscond(i,07)
      overlap(inddyy(i,5),2)=overlap(inddyy(i,5),2)+fscond(i,08)
      overlap(inddyy(i,6),2)=overlap(inddyy(i,6),2)+fscond(i,09)
      overlap(inddyy(i,1),3)=overlap(inddyy(i,1),3)+fscond(i,03)
      overlap(inddyy(i,2),3)=overlap(inddyy(i,2),3)+fscond(i,05)
      overlap(inddyy(i,3),3)=overlap(inddyy(i,3),3)+fscond(i,06)
      overlap(inddyy(i,4),3)=overlap(inddyy(i,4),3)+fscond(i,08)
      overlap(inddyy(i,5),3)=overlap(inddyy(i,5),3)+fscond(i,09)
 1013 overlap(inddyy(i,6),3)=overlap(inddyy(i,6),3)+fscond(i,10)

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

      do 1015 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)
 1015 ps(i,3)=pzaz(i)*ss(i)+eta2(i)*ps0(i,3)

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

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

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

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

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

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

      do 1022 j=1,3
      do 1022 i=1,n
 1022 pscond(icontr(i),j)=pscond(icontr(i),j)+ps(i,j)

      do 1023 j=1,10
      do 1023 i=1,n
 1023 fscond(icontr(i),j)=fscond(icontr(i),j)+fs(i,j)*twozetaa(i)

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

      do 1025 k=1,3
      do 1025 j=1,3
      do 1025 i=1,m
 1025 gradient(k,icentr(i,1))=gradient(k,icentr(i,1))-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)

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

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

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

      return
      end
