      subroutine coreps(m,n,ncentr,icontr,inddyy,iloc,coord,charge,
     &                  overlap,core,zeta,px,py,pz,pxax,pyay,pzaz,eta,
     &                  base,etarr,ss,ss0,ss1,ps,pscond,pxcx,pycy,pzcz,
     &                  u,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 icontr(*),inddyy(m,*),iloc(*)

      dimension coord(3,*),charge(*),overlap(*),core(*),zeta(*),
     &          px(*),py(*),pz(*),pxax(*),pyay(*),pzaz(*),eta(*),
     &          base(*),etarr(*),ss(*),ss0(*),ss1(*),ps(n,*),
     &          pscond(n,*),pxcx(*),pycy(*),pzcz(*),u(*),eta2(*)

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

      do 1001 i=1,n
 1001 eta2(i)=two*eta(i)

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

      do 1003 i=1,n
      ps(i,1)=pxax(i)*ss(i)
      ps(i,2)=pyay(i)*ss(i)
 1003 ps(i,3)=pzaz(i)*ss(i)

      do 1004 j=1,3
      do 1004 i=1,n
 1004 overlap(inddyy(icontr(i),j))=overlap(inddyy(icontr(i),j))+ps(i,j)

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

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

      do 1007 j=1,3
      do 1007 i=1,n
 1007 core(inddyy(icontr(i),j))=core(inddyy(icontr(i),j))+ps(i,j)

      do 1008 j=1,3
      do 1008 i=1,n
 1008 pscond(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(2,n,iloc,u,ps,ss,ss0,ss1)

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

      do 2004 i=1,n
      ps(i,1)=pxax(i)*ss0(i)-pxcx(i)*ss1(i)
      ps(i,2)=pyay(i)*ss0(i)-pycy(i)*ss1(i)
 2004 ps(i,3)=pzaz(i)*ss0(i)-pzcz(i)*ss1(i)

      do 2005 j=1,3
      do 2005 i=1,n
 2005 pscond(i,j)=pscond(i,j)+ps(i,j)

 2001 continue

      do 2006 j=1,3
      do 2006 i=1,n
 2006 core(inddyy(icontr(i),j))=core(inddyy(icontr(i),j))-pscond(i,j)

      return
      end
