      subroutine gradovpp(ndim,m,n,iminj,icontr,inddyy,icentr,density,
     &                    overlap,gradient,zeta,twozetaa,eta,eta2,base,
     &                    etarr,pxax,pyay,pzaz,pxbx,pyby,pzbz,twozeta,
     &                    ss0,ss,ps,sp,ds,dp,c,c1,spcond,dpcond,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 iminj(*),icontr(*),inddyy(m,*),icentr(m,*)

      dimension density(*),overlap(ndim,*),gradient(3,*),zeta(*),
     &          twozetaa(*),eta(*),eta2(*),base(*),etarr(*),pxax(*),
     &          pyay(*),pzaz(*),pxbx(*),pyby(*),pzbz(*),twozeta(*),
     &          ss0(*),ss(*),ps(n,*),sp(n,*),ds(n,*),dp(n,*),c(*),
     &          c1(n,*),spcond(m,*),dpcond(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
      ps(i,1)=pxax(i)*ss0(i)
      ps(i,2)=pyay(i)*ss0(i)
 1003 ps(i,3)=pzaz(i)*ss0(i)

      do 1004 i=1,n
      sp(i,1)=pxbx(i)*ss0(i)
      sp(i,2)=pyby(i)*ss0(i)
 1004 sp(i,3)=pzbz(i)*ss0(i)

      do 1005 i=1,n
 1005 c(i)=ss0(i)/twozeta(i)

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

      do 1007 j=1,3
      do 1007 i=1,n
 1007 c1(i,j)=ps(i,j)/twozeta(i)

      do 1008 i=1,n
      dp(i,01)=pxbx(i)*ds(i,1)+c1(i,1)*two
      dp(i,02)=pyby(i)*ds(i,1)
      dp(i,03)=pzbz(i)*ds(i,1)
      dp(i,04)=pxbx(i)*ds(i,2)+c1(i,2)
      dp(i,05)=pyby(i)*ds(i,2)+c1(i,1)
      dp(i,06)=pzbz(i)*ds(i,2)
      dp(i,07)=pxbx(i)*ds(i,3)+c1(i,3)
      dp(i,08)=pyby(i)*ds(i,3)
      dp(i,09)=pzbz(i)*ds(i,3)+c1(i,1)
      dp(i,10)=pxbx(i)*ds(i,4)
      dp(i,11)=pyby(i)*ds(i,4)+c1(i,2)*two
      dp(i,12)=pzbz(i)*ds(i,4)
      dp(i,13)=pxbx(i)*ds(i,5)
      dp(i,14)=pyby(i)*ds(i,5)+c1(i,3)
      dp(i,15)=pzbz(i)*ds(i,5)+c1(i,2)
      dp(i,16)=pxbx(i)*ds(i,6)
      dp(i,17)=pyby(i)*ds(i,6)
 1008 dp(i,18)=pzbz(i)*ds(i,6)+c1(i,3)*two

      do 1009 j=1,3
      do 1009 i=1,m
 1009 spcond(i,j)=zero

      do 1010 j=1,18
      do 1010 i=1,m
 1010 dpcond(i,j)=zero

      do 1011 j=1,3
      do 1011 i=1,n
 1011 spcond(icontr(i),j)=spcond(icontr(i),j)+sp(i,j)

      do 1012 j=1,18
      do 1012 i=1,n
 1012 dpcond(icontr(i),j)=dpcond(icontr(i),j)+dp(i,j)*twozetaa(i)

      do 1013 i=1,m
      overlap(inddyy(i,1),1)=overlap(inddyy(i,1),1)-spcond(i,1)
      overlap(inddyy(i,4),2)=overlap(inddyy(i,4),2)-spcond(i,1)
      overlap(inddyy(i,5),2)=overlap(inddyy(i,5),2)-spcond(i,2)
      overlap(inddyy(i,7),3)=overlap(inddyy(i,7),3)-spcond(i,1)
      overlap(inddyy(i,8),3)=overlap(inddyy(i,8),3)-spcond(i,2)
 1013 overlap(inddyy(i,9),3)=overlap(inddyy(i,9),3)-spcond(i,3)

      do 1014 i=1,m
      if(iminj(i).ne.0) then

              overlap(inddyy(i,2),1)=overlap(inddyy(i,2),1)-spcond(i,2)
              overlap(inddyy(i,3),1)=overlap(inddyy(i,3),1)-spcond(i,3)
              overlap(inddyy(i,6),2)=overlap(inddyy(i,6),2)-spcond(i,3)

                        endif
 1014 continue

      do 1015 i=1,m
      overlap(inddyy(i,1),1)=overlap(inddyy(i,1),1)+dpcond(i,01)
      overlap(inddyy(i,4),1)=overlap(inddyy(i,4),1)+dpcond(i,04)
      overlap(inddyy(i,5),1)=overlap(inddyy(i,5),1)+dpcond(i,05)
      overlap(inddyy(i,7),1)=overlap(inddyy(i,7),1)+dpcond(i,07)
      overlap(inddyy(i,8),1)=overlap(inddyy(i,8),1)+dpcond(i,08)
      overlap(inddyy(i,9),1)=overlap(inddyy(i,9),1)+dpcond(i,09)
      overlap(inddyy(i,1),2)=overlap(inddyy(i,1),2)+dpcond(i,04)
      overlap(inddyy(i,4),2)=overlap(inddyy(i,4),2)+dpcond(i,10)
      overlap(inddyy(i,5),2)=overlap(inddyy(i,5),2)+dpcond(i,11)
      overlap(inddyy(i,7),2)=overlap(inddyy(i,7),2)+dpcond(i,13)
      overlap(inddyy(i,8),2)=overlap(inddyy(i,8),2)+dpcond(i,14)
      overlap(inddyy(i,9),2)=overlap(inddyy(i,9),2)+dpcond(i,15)
      overlap(inddyy(i,1),3)=overlap(inddyy(i,1),3)+dpcond(i,07)
      overlap(inddyy(i,4),3)=overlap(inddyy(i,4),3)+dpcond(i,13)
      overlap(inddyy(i,5),3)=overlap(inddyy(i,5),3)+dpcond(i,14)
      overlap(inddyy(i,7),3)=overlap(inddyy(i,7),3)+dpcond(i,16)
      overlap(inddyy(i,8),3)=overlap(inddyy(i,8),3)+dpcond(i,17)
 1015 overlap(inddyy(i,9),3)=overlap(inddyy(i,9),3)+dpcond(i,18)

      do 1016 i=1,m
      if(iminj(i).ne.0) then

             overlap(inddyy(i,2),1)=overlap(inddyy(i,2),1)+dpcond(i,02)
             overlap(inddyy(i,3),1)=overlap(inddyy(i,3),1)+dpcond(i,03)
             overlap(inddyy(i,6),1)=overlap(inddyy(i,6),1)+dpcond(i,06)
             overlap(inddyy(i,2),2)=overlap(inddyy(i,2),2)+dpcond(i,05)
             overlap(inddyy(i,3),2)=overlap(inddyy(i,3),2)+dpcond(i,06)
             overlap(inddyy(i,6),2)=overlap(inddyy(i,6),2)+dpcond(i,12)
             overlap(inddyy(i,2),3)=overlap(inddyy(i,2),3)+dpcond(i,08)
             overlap(inddyy(i,3),3)=overlap(inddyy(i,3),3)+dpcond(i,09)
             overlap(inddyy(i,6),3)=overlap(inddyy(i,6),3)+dpcond(i,15)

                        endif
 1016 continue

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

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

      do 1019 i=1,n
      sp(i,1)=pxbx(i)*ss(i)+eta2(i)*sp(i,1)
      sp(i,2)=pyby(i)*ss(i)+eta2(i)*sp(i,2)
 1019 sp(i,3)=pzbz(i)*ss(i)+eta2(i)*sp(i,3)

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

      do 1021 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)
 1021 ds(i,6)=pzaz(i)*ps(i,3)+eta2(i)*ds(i,6)+c(i)

      do 1022 j=1,3
      do 1022 i=1,n
 1022 c1(i,j)=ps(i,j)/twozeta(i)

      do 1023 i=1,n
      dp(i,01)=pxbx(i)*ds(i,1)+eta2(i)*dp(i,01)+c1(i,1)*two
      dp(i,02)=pyby(i)*ds(i,1)+eta2(i)*dp(i,02)
      dp(i,03)=pzbz(i)*ds(i,1)+eta2(i)*dp(i,03)
      dp(i,04)=pxbx(i)*ds(i,2)+eta2(i)*dp(i,04)+c1(i,2)
      dp(i,05)=pyby(i)*ds(i,2)+eta2(i)*dp(i,05)+c1(i,1)
      dp(i,06)=pzbz(i)*ds(i,2)+eta2(i)*dp(i,06)
      dp(i,07)=pxbx(i)*ds(i,3)+eta2(i)*dp(i,07)+c1(i,3)
      dp(i,08)=pyby(i)*ds(i,3)+eta2(i)*dp(i,08)
      dp(i,09)=pzbz(i)*ds(i,3)+eta2(i)*dp(i,09)+c1(i,1)
      dp(i,10)=pxbx(i)*ds(i,4)+eta2(i)*dp(i,10)
      dp(i,11)=pyby(i)*ds(i,4)+eta2(i)*dp(i,11)+c1(i,2)*two
      dp(i,12)=pzbz(i)*ds(i,4)+eta2(i)*dp(i,12)
      dp(i,13)=pxbx(i)*ds(i,5)+eta2(i)*dp(i,13)
      dp(i,14)=pyby(i)*ds(i,5)+eta2(i)*dp(i,14)+c1(i,3)
      dp(i,15)=pzbz(i)*ds(i,5)+eta2(i)*dp(i,15)+c1(i,2)
      dp(i,16)=pxbx(i)*ds(i,6)+eta2(i)*dp(i,16)
      dp(i,17)=pyby(i)*ds(i,6)+eta2(i)*dp(i,17)
 1023 dp(i,18)=pzbz(i)*ds(i,6)+eta2(i)*dp(i,18)+c1(i,3)*two

      do 1024 j=1,3
      do 1024 i=1,m
 1024 spcond(i,j)=zero

      do 1025 j=1,18
      do 1025 i=1,m
 1025 dpcond(i,j)=zero

      do 1026 j=1,3
      do 1026 i=1,n
 1026 spcond(icontr(i),j)=spcond(icontr(i),j)+sp(i,j)

      do 1027 j=1,18
      do 1027 i=1,n
 1027 dpcond(icontr(i),j)=dpcond(icontr(i),j)+dp(i,j)*twozetaa(i)

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

      do 1029 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,2)=zero
                              grad(i,3)=zero
                              grad(i,6)=zero
                        endif
 1029 continue

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

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

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

      do 1033 i=1,m
      if(iminj(i).eq.0) then
                              grad(i,02)=zero
                              grad(i,03)=zero
                              grad(i,06)=zero
                              grad(i,11)=zero
                              grad(i,12)=zero
                              grad(i,15)=zero
                              grad(i,20)=zero
                              grad(i,21)=zero
                              grad(i,24)=zero
                        endif
 1033 continue

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

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

      return
      end
