      subroutine boys_ps(m,n,icontr,inddyy,zeta,px,py,pz,pxax,pyay,pzaz,
     &                   base,ss,ps,c,x_moment,y_moment,z_moment)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      dimension icontr(*),inddyy(m,*)

      dimension zeta(*),px(*),py(*),pz(*),pxax(*),pyay(*),pzaz(*),
     &          base(*),ss(n,*),ps(n,*),c(*),x_moment(*),y_moment(*),
     &          z_moment(*)

      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
      ss(i,1)=px(i)*base(i)
      ss(i,2)=py(i)*base(i)
 1002 ss(i,3)=pz(i)*base(i)

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

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

      do 1005 j=1,3
      do 1005 i=1,n
      ii=inddyy(icontr(i),j)
      x_moment(ii)=x_moment(ii)+ps(i,3*(j-1)+1)
      y_moment(ii)=y_moment(ii)+ps(i,3*(j-1)+2)
 1005 z_moment(ii)=z_moment(ii)+ps(i,3*(j-1)+3)

      return
      end
