      subroutine dipoleds(m,n,icontr,inddyy,density,zeta,px,py,pz,pxax,
     &                    pyay,pzaz,base,ss,ps,ds,pso,c,c1,c2,dscond,
     &                    dipole)

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,*)

      dimension density(*),zeta(*),px(*),py(*),pz(*),pxax(*),pyay(*),
     &          pzaz(*),base(*),ss(n,*),ps(n,*),ds(n,*),pso(n,*),c(*),
     &          c1(n,*),c2(n,*),dscond(m,*),dipole(*)

      data zero,thrhlf,two,pi/0.0d0,1.5d0,2.0d0,3.1415926535898d0/

      do 1001 i=1,n
 1001 base(i)=base(i)*((pi/zeta(i))**thrhlf)

      do 1002 i=1,n
      pso(i,1)=pxax(i)*base(i)
      pso(i,2)=pyay(i)*base(i)
 1002 pso(i,3)=pzaz(i)*base(i)

      do 1003 i=1,n
      ss(i,1)=px(i)*base(i)
      ss(i,2)=py(i)*base(i)
 1003 ss(i,3)=pz(i)*base(i)

      do 1004 i=1,n
 1004 c(i)=base(i)/(two*zeta(i))

      do 1005 i=1,n
      ps(i,1)=pxax(i)*ss(i,1)+c(i)
      ps(i,2)=pxax(i)*ss(i,2)
      ps(i,3)=pxax(i)*ss(i,3)
      ps(i,4)=pyay(i)*ss(i,1)
      ps(i,5)=pyay(i)*ss(i,2)+c(i)
      ps(i,6)=pyay(i)*ss(i,3)
      ps(i,7)=pzaz(i)*ss(i,1)
      ps(i,8)=pzaz(i)*ss(i,2)
 1005 ps(i,9)=pzaz(i)*ss(i,3)+c(i)

      do 1006 j=1,3
      do 1006 i=1,n
 1006 c1(i,j)=ss(i,j)/(two*zeta(i))

      do 1007 j=1,3
      do 1007 i=1,n
 1007 c2(i,j)=pso(i,j)/(two*zeta(i))

      do 1008 i=1,n
      ds(i,01)=pxax(i)*ps(i,1)+c1(i,1)+c2(i,1)
      ds(i,02)=pxax(i)*ps(i,2)+c1(i,2)
      ds(i,03)=pxax(i)*ps(i,3)+c1(i,3)
      ds(i,04)=pyay(i)*ps(i,1)
      ds(i,05)=pyay(i)*ps(i,2)        +c2(i,1)
      ds(i,06)=pyay(i)*ps(i,3)
      ds(i,07)=pzaz(i)*ps(i,1)
      ds(i,08)=pzaz(i)*ps(i,2)
      ds(i,09)=pzaz(i)*ps(i,3)        +c2(i,1)
      ds(i,10)=pyay(i)*ps(i,4)+c1(i,1)
      ds(i,11)=pyay(i)*ps(i,5)+c1(i,2)+c2(i,2)
      ds(i,12)=pyay(i)*ps(i,6)+c1(i,3)
      ds(i,13)=pzaz(i)*ps(i,4)
      ds(i,14)=pzaz(i)*ps(i,5)
      ds(i,15)=pzaz(i)*ps(i,6)        +c2(i,2)
      ds(i,16)=pzaz(i)*ps(i,7)+c1(i,1)
      ds(i,17)=pzaz(i)*ps(i,8)+c1(i,2)
 1008 ds(i,18)=pzaz(i)*ps(i,9)+c1(i,3)+c2(i,3)

      do 1009 j=1,18
      do 1009 i=1,m
 1009 dscond(i,j)=zero

      do 1010 j=1,18
      do 1010 i=1,n
 1010 dscond(icontr(i),j)=dscond(icontr(i),j)+ds(i,j)

      do 1011 k=1,6
      do 1011 j=1,3
      do 1011 i=1,m
 1011 dipole(j)=dipole(j)-density(inddyy(i,k))*dscond(i,3*(k-1)+j)

      return
      end
