      subroutine corepp(m,n,ncentr,iminj,icontr,inddyy,iloc,coord,
     &                  charge,overlap,core,zeta,axbx,ayby,azbz,px,py,
     &                  pz,pxax,pyay,pzaz,pxbx,pyby,pzbz,eta,base,
     &                  etarr,ss,ss0,ss1,ss2,ps,ps0,ps1,ds0,pp,ppcond,
     &                  c,pxcx,pycy,pzcz,u,twozeta,eta2)


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

      dimension coord(3,*),charge(*),overlap(*),core(*),zeta(*),
     &          axbx(*),ayby(*),azbz(*),px(*),py(*),pz(*),pxax(*),
     &          pyay(*),pzaz(*),pxbx(*),pyby(*),pzbz(*),eta(*),base(*),
     &          etarr(*),ss(*),ss0(*),ss1(*),ss2(*),ps(n,*),ps0(n,*),
     &          ps1(n,*),ds0(n,*),pp(n,*),ppcond(n,*),c(*),pxcx(*),
     &          pycy(*),pzcz(*),u(*),twozeta(*),eta2(*)

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

      do 1001 i=1,n
      eta2(i)=two*eta(i)
 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
      pp(i,1)=pxbx(i)*ps0(i,1)+c(i)
      pp(i,2)=pyby(i)*ps0(i,1)
      pp(i,3)=pzbz(i)*ps0(i,1)
      pp(i,4)=pxbx(i)*ps0(i,2)
      pp(i,5)=pyby(i)*ps0(i,2)+c(i)
      pp(i,6)=pzbz(i)*ps0(i,2)
      pp(i,7)=pxbx(i)*ps0(i,3)
      pp(i,8)=pyby(i)*ps0(i,3)
 1005 pp(i,9)=pzbz(i)*ps0(i,3)+c(i)

      do 1006 i=1,n
      overlap(inddyy(icontr(i),1))=overlap(inddyy(icontr(i),1))+pp(i,1)
      overlap(inddyy(icontr(i),4))=overlap(inddyy(icontr(i),4))+pp(i,4)
      overlap(inddyy(icontr(i),5))=overlap(inddyy(icontr(i),5))+pp(i,5)
      overlap(inddyy(icontr(i),7))=overlap(inddyy(icontr(i),7))+pp(i,7)
      overlap(inddyy(icontr(i),8))=overlap(inddyy(icontr(i),8))+pp(i,8)
 1006 overlap(inddyy(icontr(i),9))=overlap(inddyy(icontr(i),9))+pp(i,9)


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

              overlap(inddyy(icontr(i),2))=overlap(inddyy(icontr(i),2))
     &                                    +pp(i,2)
              overlap(inddyy(icontr(i),3))=overlap(inddyy(icontr(i),3))
     &                                    +pp(i,3)
              overlap(inddyy(icontr(i),6))=overlap(inddyy(icontr(i),6))
     &                                    +pp(i,6)

                        endif
 1007 continue

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

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

      do 1010 i=1,n
 1010 c(i)=ss(i)/twozeta(i)

      do 1011 i=1,n
      pp(i,1)=pxbx(i)*ps(i,1)+eta2(i)*pp(i,1)+c(i)
      pp(i,2)=pyby(i)*ps(i,1)+eta2(i)*pp(i,2)
      pp(i,3)=pzbz(i)*ps(i,1)+eta2(i)*pp(i,3)
      pp(i,4)=pxbx(i)*ps(i,2)+eta2(i)*pp(i,4)
      pp(i,5)=pyby(i)*ps(i,2)+eta2(i)*pp(i,5)+c(i)
      pp(i,6)=pzbz(i)*ps(i,2)+eta2(i)*pp(i,6)
      pp(i,7)=pxbx(i)*ps(i,3)+eta2(i)*pp(i,7)
      pp(i,8)=pyby(i)*ps(i,3)+eta2(i)*pp(i,8)
 1011 pp(i,9)=pzbz(i)*ps(i,3)+eta2(i)*pp(i,9)+c(i)

      do 1012 i=1,n
      core(inddyy(icontr(i),1))=core(inddyy(icontr(i),1))+pp(i,1)
      core(inddyy(icontr(i),4))=core(inddyy(icontr(i),4))+pp(i,4)
      core(inddyy(icontr(i),5))=core(inddyy(icontr(i),5))+pp(i,5)
      core(inddyy(icontr(i),7))=core(inddyy(icontr(i),7))+pp(i,7)
      core(inddyy(icontr(i),8))=core(inddyy(icontr(i),8))+pp(i,8)
 1012 core(inddyy(icontr(i),9))=core(inddyy(icontr(i),9))+pp(i,9)

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

               core(inddyy(icontr(i),2))=core(inddyy(icontr(i),2))
     &                                  +pp(i,2)
               core(inddyy(icontr(i),3))=core(inddyy(icontr(i),3))
     &                                  +pp(i,3)
               core(inddyy(icontr(i),6))=core(inddyy(icontr(i),6))
     &                                  +pp(i,6)

                        endif
 1013 continue

      do 1014 j=1,9
      do 1014 i=1,n
 1014 ppcond(i,j)=zero

      do 2001 nn=1,ncentr

      ccx=coord(1,nn)
      ccy=coord(2,nn)
      ccz=coord(3,nn)

      twopichrg=two*pi*charge(nn)

      do 2002 i=1,n
      pxcx(i)=px(i)-ccx
      pycy(i)=py(i)-ccy
      pzcz(i)=pz(i)-ccz
 2002 u(i)=zeta(i)*(pxcx(i)*pxcx(i)+pycy(i)*pycy(i)+pzcz(i)*pzcz(i))

      call auggv(3,n,iloc,u,ps0,ss0,eta2,c)

      do 2003 i=1,n
      ss0(i)=base(i)*ps0(i,1)*twopichrg/zeta(i)
      ss1(i)=base(i)*ps0(i,2)*twopichrg/zeta(i)
 2003 ss2(i)=base(i)*ps0(i,3)*twopichrg/zeta(i)

      do 2004 i=1,n

      ps0(i,1)=pxax(i)*ss0(i)-pxcx(i)*ss1(i)
      ps0(i,2)=pyay(i)*ss0(i)-pycy(i)*ss1(i)
      ps0(i,3)=pzaz(i)*ss0(i)-pzcz(i)*ss1(i)

      ps1(i,1)=pxax(i)*ss1(i)-pxcx(i)*ss2(i)
      ps1(i,2)=pyay(i)*ss1(i)-pycy(i)*ss2(i)
      ps1(i,3)=pzaz(i)*ss1(i)-pzcz(i)*ss2(i)

 2004 continue

      do 2005 i=1,n
 2005 c(i)=(ss0(i)-ss1(i))/twozeta(i)

      do 2006 i=1,n
      ds0(i,1)=pxax(i)*ps0(i,1)-pxcx(i)*ps1(i,1)+c(i)
      ds0(i,2)=pyay(i)*ps0(i,1)-pycy(i)*ps1(i,1)
      ds0(i,3)=pzaz(i)*ps0(i,1)-pzcz(i)*ps1(i,1)
      ds0(i,4)=pyay(i)*ps0(i,2)-pycy(i)*ps1(i,2)+c(i)
      ds0(i,5)=pzaz(i)*ps0(i,2)-pzcz(i)*ps1(i,2)
 2006 ds0(i,6)=pzaz(i)*ps0(i,3)-pzcz(i)*ps1(i,3)+c(i)

      do 2011 i=1,n
      ppcond(i,1)=ppcond(i,1)+axbx(icontr(i))*ps0(i,1)+ds0(i,1)
      ppcond(i,2)=ppcond(i,2)+ayby(icontr(i))*ps0(i,1)+ds0(i,2)
      ppcond(i,3)=ppcond(i,3)+azbz(icontr(i))*ps0(i,1)+ds0(i,3)
      ppcond(i,4)=ppcond(i,4)+axbx(icontr(i))*ps0(i,2)+ds0(i,2)
      ppcond(i,5)=ppcond(i,5)+ayby(icontr(i))*ps0(i,2)+ds0(i,4)
      ppcond(i,6)=ppcond(i,6)+azbz(icontr(i))*ps0(i,2)+ds0(i,5)
      ppcond(i,7)=ppcond(i,7)+axbx(icontr(i))*ps0(i,3)+ds0(i,3)
      ppcond(i,8)=ppcond(i,8)+ayby(icontr(i))*ps0(i,3)+ds0(i,5)
 2011 ppcond(i,9)=ppcond(i,9)+azbz(icontr(i))*ps0(i,3)+ds0(i,6)

 2001 continue

      do 2012 i=1,n
      core(inddyy(icontr(i),1))=core(inddyy(icontr(i),1))-ppcond(i,1)
      core(inddyy(icontr(i),4))=core(inddyy(icontr(i),4))-ppcond(i,4)
      core(inddyy(icontr(i),5))=core(inddyy(icontr(i),5))-ppcond(i,5)
      core(inddyy(icontr(i),7))=core(inddyy(icontr(i),7))-ppcond(i,7)
      core(inddyy(icontr(i),8))=core(inddyy(icontr(i),8))-ppcond(i,8)
 2012 core(inddyy(icontr(i),9))=core(inddyy(icontr(i),9))-ppcond(i,9)

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

               core(inddyy(icontr(i),2))=core(inddyy(icontr(i),2))
     &                                  -ppcond(i,2)
               core(inddyy(icontr(i),3))=core(inddyy(icontr(i),3))
     &                                  -ppcond(i,3)
               core(inddyy(icontr(i),6))=core(inddyy(icontr(i),6))
     &                                  -ppcond(i,6)

                        endif
 2013 continue

      return
      end
