      subroutine coreds(m,n,ncentr,icontr,inddyy,iloc,coord,charge,
     &                  overlap,core,zeta,zetaa,px,py,pz,pxax,pyay,pzaz,
     &                  eta,base,etarr,ss,ss0,ss1,ss2,ps,ps0,ps1,ds,
     &                  dscond,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 icontr(*),inddyy(m,*),iloc(*)

      dimension coord(3,*),charge(*),overlap(*),core(*),zeta(*),
     &          zetaa(*),px(*),py(*),pz(*),pxax(*),pyay(*),pzaz(*),
     &          eta(*),base(*),etarr(*),ss(*),ss0(*),ss1(*),ss2(*),
     &          ps(n,*),ps0(n,*),ps1(n,*),ds(n,*),dscond(n,*),c(*),
     &          pxcx(*),pycy(*),pzcz(*),u(*),twozeta(*),eta2(*)

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

      sqrt3=sqrt(three)

      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
      ds(i,1)=pxax(i)*ps0(i,1)+c(i)
      ds(i,2)=pyay(i)*ps0(i,1)
      ds(i,3)=pzaz(i)*ps0(i,1)
      ds(i,4)=pyay(i)*ps0(i,2)+c(i)
      ds(i,5)=pzaz(i)*ps0(i,2)
 1005 ds(i,6)=pzaz(i)*ps0(i,3)+c(i)

      do 1006 j=1,6
      do 1006 i=1,n
 1006 overlap(inddyy(icontr(i),j))=overlap(inddyy(icontr(i),j))+ds(i,j)

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

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

      do 1009 i=1,n
 1009 c(i)=ss(i)/twozeta(i)-eta2(i)*ss0(i)/(two*zetaa(i))

      do 1010 i=1,n
      ds(i,1)=pxax(i)*ps(i,1)+eta2(i)*ds(i,1)+c(i)
      ds(i,2)=pyay(i)*ps(i,1)+eta2(i)*ds(i,2)
      ds(i,3)=pzaz(i)*ps(i,1)+eta2(i)*ds(i,3)
      ds(i,4)=pyay(i)*ps(i,2)+eta2(i)*ds(i,4)+c(i)
      ds(i,5)=pzaz(i)*ps(i,2)+eta2(i)*ds(i,5)
 1010 ds(i,6)=pzaz(i)*ps(i,3)+eta2(i)*ds(i,6)+c(i)

      do 1011 j=1,6
      do 1011 i=1,n
 1011 core(inddyy(icontr(i),j))=core(inddyy(icontr(i),j))+ds(i,j)

      do 1012 j=1,6
      do 1012 i=1,n
 1012 dscond(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,ss1,ss2)

      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
      ds(i,1)=pxax(i)*ps0(i,1)-pxcx(i)*ps1(i,1)+c(i)
      ds(i,2)=pyay(i)*ps0(i,1)-pycy(i)*ps1(i,1)
      ds(i,3)=pzaz(i)*ps0(i,1)-pzcz(i)*ps1(i,1)
      ds(i,4)=pyay(i)*ps0(i,2)-pycy(i)*ps1(i,2)+c(i)
      ds(i,5)=pzaz(i)*ps0(i,2)-pzcz(i)*ps1(i,2)
 2006 ds(i,6)=pzaz(i)*ps0(i,3)-pzcz(i)*ps1(i,3)+c(i)

      do 2007 j=1,6
      do 2007 i=1,n
 2007 dscond(i,j)=dscond(i,j)+ds(i,j)

 2001 continue

      do 2008 j=1,6
      do 2008 i=1,n
 2008 core(inddyy(icontr(i),j))=core(inddyy(icontr(i),j))-dscond(i,j)

      return
      end
