      subroutine boys_pp(m,n,iminj,icontr,inddyy,zeta,px,py,pz,pxax,
     &                   pyay,pzaz,pxbx,pyby,pzbz,base,ss,ps,pp,pso,c,
     &                   c1,c2,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 iminj(*),icontr(*),inddyy(m,*)

      dimension zeta(*),px(*),py(*),pz(*),pxax(*),pyay(*),pzaz(*),
     &          pxbx(*),pyby(*),pzbz(*),base(*),ss(n,*),ps(n,*),pp(n,*),
     &          pso(n,*),c(*),c1(n,*),c2(n,*),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
      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
      pp(i,01)=pxbx(i)*ps(i,1)+c1(i,1)+c2(i,1)
      pp(i,02)=pxbx(i)*ps(i,2)+c1(i,2)
      pp(i,03)=pxbx(i)*ps(i,3)+c1(i,3)
      pp(i,04)=pyby(i)*ps(i,1)
      pp(i,05)=pyby(i)*ps(i,2)        +c2(i,1)
      pp(i,06)=pyby(i)*ps(i,3)
      pp(i,07)=pzbz(i)*ps(i,1)
      pp(i,08)=pzbz(i)*ps(i,2)
      pp(i,09)=pzbz(i)*ps(i,3)        +c2(i,1)
      pp(i,10)=pxbx(i)*ps(i,4)        +c2(i,2)
      pp(i,11)=pxbx(i)*ps(i,5)
      pp(i,12)=pxbx(i)*ps(i,6)
      pp(i,13)=pyby(i)*ps(i,4)+c1(i,1)
      pp(i,14)=pyby(i)*ps(i,5)+c1(i,2)+c2(i,2)
      pp(i,15)=pyby(i)*ps(i,6)+c1(i,3)
      pp(i,16)=pzbz(i)*ps(i,4)
      pp(i,17)=pzbz(i)*ps(i,5)
      pp(i,18)=pzbz(i)*ps(i,6)        +c2(i,2)
      pp(i,19)=pxbx(i)*ps(i,7)        +c2(i,3)
      pp(i,20)=pxbx(i)*ps(i,8)
      pp(i,21)=pxbx(i)*ps(i,9)
      pp(i,22)=pyby(i)*ps(i,7)
      pp(i,23)=pyby(i)*ps(i,8)        +c2(i,3)
      pp(i,24)=pyby(i)*ps(i,9)
      pp(i,25)=pzbz(i)*ps(i,7)+c1(i,1)
      pp(i,26)=pzbz(i)*ps(i,8)+c1(i,2)
 1008 pp(i,27)=pzbz(i)*ps(i,9)+c1(i,3)+c2(i,3)

      do 1009 i=1,n
      ii=icontr(i)
      x_moment(inddyy(ii,1))=x_moment(inddyy(ii,1))+pp(i,01)
      x_moment(inddyy(ii,4))=x_moment(inddyy(ii,4))+pp(i,10)
      x_moment(inddyy(ii,5))=x_moment(inddyy(ii,5))+pp(i,13)
      x_moment(inddyy(ii,7))=x_moment(inddyy(ii,7))+pp(i,19)
      x_moment(inddyy(ii,8))=x_moment(inddyy(ii,8))+pp(i,22)
      x_moment(inddyy(ii,9))=x_moment(inddyy(ii,9))+pp(i,25)
      y_moment(inddyy(ii,1))=y_moment(inddyy(ii,1))+pp(i,02)
      y_moment(inddyy(ii,4))=y_moment(inddyy(ii,4))+pp(i,11)
      y_moment(inddyy(ii,5))=y_moment(inddyy(ii,5))+pp(i,14)
      y_moment(inddyy(ii,7))=y_moment(inddyy(ii,7))+pp(i,20)
      y_moment(inddyy(ii,8))=y_moment(inddyy(ii,8))+pp(i,23)
      y_moment(inddyy(ii,9))=y_moment(inddyy(ii,9))+pp(i,26)
      z_moment(inddyy(ii,1))=z_moment(inddyy(ii,1))+pp(i,03)
      z_moment(inddyy(ii,4))=z_moment(inddyy(ii,4))+pp(i,12)
      z_moment(inddyy(ii,5))=z_moment(inddyy(ii,5))+pp(i,15)
      z_moment(inddyy(ii,7))=z_moment(inddyy(ii,7))+pp(i,21)
      z_moment(inddyy(ii,8))=z_moment(inddyy(ii,8))+pp(i,24)
 1009 z_moment(inddyy(ii,9))=z_moment(inddyy(ii,9))+pp(i,27)

      do 1010 i=1,n
      ii=icontr(i)
      if(iminj(ii).ne.0) then

               x_moment(inddyy(ii,2))=x_moment(inddyy(ii,2))+pp(i,04)
               x_moment(inddyy(ii,3))=x_moment(inddyy(ii,3))+pp(i,07)
               x_moment(inddyy(ii,6))=x_moment(inddyy(ii,6))+pp(i,16)
               y_moment(inddyy(ii,2))=y_moment(inddyy(ii,2))+pp(i,05)
               y_moment(inddyy(ii,3))=y_moment(inddyy(ii,3))+pp(i,08)
               y_moment(inddyy(ii,6))=y_moment(inddyy(ii,6))+pp(i,17)
               z_moment(inddyy(ii,2))=z_moment(inddyy(ii,2))+pp(i,06)
               z_moment(inddyy(ii,3))=z_moment(inddyy(ii,3))+pp(i,09)
               z_moment(inddyy(ii,6))=z_moment(inddyy(ii,6))+pp(i,18)

                         endif
 1010 continue

      return
      end
