      subroutine chrgss(m,n,ncds,ncdspd,ilfunc,icfunc,icontr,iicontr,
     &                  inddyy,iloc,g,coord,cdfitc,alphaaux,coeffs,
     &                  coeffp,coeffd,fock,zeta,px,py,pz,base,u,zetarec,
     &                  wxcx,wycy,wzcz,sss,sss0,sss1,sss2,ssp,ssp0,ssp1,
     &                  ssd,ssd0,zetac,c,zzeta,ppx,ppy,ppz,gather)

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 cdfitp(3),cdfitd(6)

      dimension ilfunc(*),icfunc(*),icontr(*),iicontr(*),inddyy(*),
     &          iloc(*)

      dimension g(*),coord(3,*),cdfitc(*),alphaaux(*),coeffs(*),
     &          coeffp(*),coeffd(*),fock(*),zeta(*),px(*),py(*),pz(*),
     &          base(*),u(*),zetarec(*),wxcx(*),wycy(*),wzcz(*),sss(*),
     &          sss0(*),sss1(*),sss2(*),ssp(m,*),ssp0(n,*),ssp1(n,*),
     &          ssd(m,*),ssd0(n,*),zetac(*),c(*),zzeta(*),ppx(*),ppy(*),
     &          ppz(*),gather(*)

      data one,two,three/1.0d0,2.0d0,3.0d0/
      data four,five,seven/4.0d0,5.0d0,7.0d0/
      data zero,half,thrhlf,pi/0.0d0,0.5d0,1.5d0,3.1415926535898d0/
      data calcut/1.0d-16/

      twoecnst=two*pi**(five/two)

      sqrt3=sqrt(three)

      dnormcd=seven/sqrt3

      do 1001 i=1,n
 1001 gather(i)=zero

      do 1002 naux=1,ncds

      ilaux=ilfunc(naux)

      ccx=coord(1,icfunc(ilaux))
      ccy=coord(2,icfunc(ilaux))
      ccz=coord(3,icfunc(ilaux))

      zetaux=alphaaux(ilaux)

      cdfits=cdfitc(naux)*coeffs(ilaux)

      do 1003 i=1,n
 1003 zetarec(i)=one/(zeta(i)+zetaux)

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

      call auggv(1,n,iloc,u,sss0,ppx,ppy,ppz)

      do 1005 i=1,n
 1005 sss0(i)=base(i)*twoecnst*sqrt(zetarec(i))*sss0(i)/zeta(i)/zetaux

      do 1006 i=1,n
 1006 gather(i)=gather(i)+sss0(i)*cdfits

 1002 continue

      do 2001 naux=1,ncdspd

      ilaux=ilfunc(ncds+naux)

      ccx=coord(1,icfunc(ilaux))
      ccy=coord(2,icfunc(ilaux))
      ccz=coord(3,icfunc(ilaux))

      zetaux=alphaaux(ilaux)

      cdfits=cdfitc(ncds+(naux-1)*10+1)*coeffs(ilaux)

      cdfitp(1)=cdfitc(ncds+(naux-1)*10+2)*coeffp(ilaux)
      cdfitp(2)=cdfitc(ncds+(naux-1)*10+3)*coeffp(ilaux)
      cdfitp(3)=cdfitc(ncds+(naux-1)*10+4)*coeffp(ilaux)

      cdfitd(1)=cdfitc(ncds+(naux-1)*10+05)*coeffd(ilaux)/dnormcd
      cdfitd(2)=cdfitc(ncds+(naux-1)*10+06)*coeffd(ilaux)
      cdfitd(3)=cdfitc(ncds+(naux-1)*10+07)*coeffd(ilaux)
      cdfitd(4)=cdfitc(ncds+(naux-1)*10+08)*coeffd(ilaux)/dnormcd
      cdfitd(5)=cdfitc(ncds+(naux-1)*10+09)*coeffd(ilaux)
      cdfitd(6)=cdfitc(ncds+(naux-1)*10+10)*coeffd(ilaux)/dnormcd

      do 2002 i=1,n
 2002 zetarec(i)=one/(zeta(i)+zetaux)

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

      call auggv(3,n,iloc,u,ssp0,ppx,ppy,ppz)

      do 2004 i=1,n
      sss0(i)=base(i)*twoecnst*sqrt(zetarec(i))*ssp0(i,1)/zeta(i)/zetaux
      sss1(i)=base(i)*twoecnst*sqrt(zetarec(i))*ssp0(i,2)/zeta(i)/zetaux
 2004 sss2(i)=base(i)*twoecnst*sqrt(zetarec(i))*ssp0(i,3)/zeta(i)/zetaux

      nn=0

      do 2005 i=1,n

      if(abs(sss0(i)).gt.calcut) then
                                       nn=nn+1

                                       iloc(nn)=i

                                       iicontr(nn)=icontr(i)

                                       zzeta(nn)=zeta(i)

                                       ppx(nn)=px(i)
                                       ppy(nn)=py(i)
                                       ppz(nn)=pz(i)

                                       zetarec(nn)=zetarec(i)

                                       sss0(nn)=sss0(i)
                                       sss1(nn)=sss1(i)
                                       sss2(nn)=sss2(i)
                                 endif

 2005 continue

      do 2006 i=1,nn
 2006 zetac(i)=zzeta(i)*zetarec(i)

      do 2007 i=1,nn
      wxcx(i)=((zzeta(i)*ppx(i)+zetaux*ccx)*zetarec(i))-ccx
      wycy(i)=((zzeta(i)*ppy(i)+zetaux*ccy)*zetarec(i))-ccy
 2007 wzcz(i)=((zzeta(i)*ppz(i)+zetaux*ccz)*zetarec(i))-ccz

      do 2008 i=1,nn

      ssp0(i,1)=wxcx(i)*sss1(i)
      ssp0(i,2)=wycy(i)*sss1(i)
      ssp0(i,3)=wzcz(i)*sss1(i)

      ssp1(i,1)=wxcx(i)*sss2(i)
      ssp1(i,2)=wycy(i)*sss2(i)
      ssp1(i,3)=wzcz(i)*sss2(i)

 2008 continue

      do 2009 i=1,nn

      c(i)=(sss0(i)-zetac(i)*sss1(i))/(two*zetaux)

      ssd0(i,1)=wxcx(i)*ssp1(i,1)+c(i)
      ssd0(i,2)=wycy(i)*ssp1(i,1)
      ssd0(i,3)=wzcz(i)*ssp1(i,1)
      ssd0(i,4)=wycy(i)*ssp1(i,2)+c(i)
      ssd0(i,5)=wzcz(i)*ssp1(i,2)
      ssd0(i,6)=wzcz(i)*ssp1(i,3)+c(i)

 2009 continue

cfpp$ nodepchk
      do 2010 i=1,nn
 2010 gather(iloc(i))=gather(iloc(i))
     &               +sss0(i)*cdfits
     &               +ssp0(i,1)*cdfitp(1)
     &               +ssp0(i,2)*cdfitp(2)
     &               +ssp0(i,3)*cdfitp(3)
     &               +ssd0(i,1)*cdfitd(1)
     &               +ssd0(i,2)*cdfitd(2)
     &               +ssd0(i,3)*cdfitd(3)
     &               +ssd0(i,4)*cdfitd(4)
     &               +ssd0(i,5)*cdfitd(5)
     &               +ssd0(i,6)*cdfitd(6)

 2001 continue

      do 1007 i=1,n
 1007 fock(inddyy(icontr(i)))=fock(inddyy(icontr(i)))+gather(i)

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

      return
      end
