      subroutine xcorss(m,n,nxcs,nxcspd,ilfunc,icfunc,icontr,iicontr,
     &                  inddyy,iloc,coord,xcfit1,xcfit2,xcfit3,alphaaux,
     &                  coeffs,coeffp,coeffd,fock1,fock2,fock3,zeta,ax,
     &                  ay,az,px,py,pz,base,zetarec,zetafac,efac,gxax,
     &                  gyay,gzaz,sss,ssscond,pss,psscond,dss,dsscond,
     &                  ssp,psp,ssd,c,zzeta,aax,aay,aaz,ppx,ppy,ppz,
     &                  axcx,aycy,azcz,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 xcfitp(3),xcfitd(6)

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

      dimension coord(3,*),xcfit1(*),xcfit2(*),xcfit3(*),alphaaux(*),
     &          coeffs(*),coeffp(*),coeffd(*),fock1(*),fock2(*),
     &          fock3(*),zeta(*),ax(*),ay(*),az(*),px(*),py(*),pz(*),
     &          base(*),zetarec(*),zetafac(*),efac(*),gxax(*),gyay(*),
     &          gzaz(*),sss(*),ssscond(*),pss(n,*),psscond(n,*),
     &          dss(n,*),dsscond(n,*),ssp(n,*),psp(n,*),ssd(n,*),c(*),
     &          zzeta(*),aax(*),aay(*),aaz(*),ppx(*),ppy(*),ppz(*),
     &          axcx(*),aycy(*),azcz(*),gather(n,*)

      data zero,one,two,three,four,five/0.0d0,1.0d0,2.0d0,
     &                                  3.0d0,4.0d0,5.0d0/
      data expcut,calcut/45.0d0,1.0d-16/

      sqrt3=sqrt(three)

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

      do 1002 naux=1,nxcs

      ilaux=ilfunc(naux)

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

      zetaux=alphaaux(ilaux)

      xcfits=xcfit3(naux)*coeffs(ilaux)

      do 1003 i=1,n
      zetafac(i)=(zeta(i)/(zeta(i)+zetaux))**(three/two)
 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 efac(i)=zetaux*zeta(i)*zetarec(i)*(pxcx*pxcx+pycy*pycy+pzcz*pzcz)

      do 1005 i=1,n
      if(efac(i).lt.expcut) then

                            sss(i)=zetafac(i)*base(i)*exp(-efac(i))
                            gather(i,1)=gather(i,1)+sss(i)*xcfits

                            endif
 1005 continue

 1002 continue

      do 2001 naux=1,nxcspd

      ilaux=ilfunc(nxcs+naux)

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

      zetaux=alphaaux(ilaux)

      xcfits=xcfit3(nxcs+(naux-1)*10+1)*coeffs(ilaux)

      xcfitp(1)=xcfit3(nxcs+(naux-1)*10+2)*coeffp(ilaux)
      xcfitp(2)=xcfit3(nxcs+(naux-1)*10+3)*coeffp(ilaux)
      xcfitp(3)=xcfit3(nxcs+(naux-1)*10+4)*coeffp(ilaux)

      xcfitd(1)=xcfit3(nxcs+(naux-1)*10+05)*coeffd(ilaux)/sqrt3
      xcfitd(2)=xcfit3(nxcs+(naux-1)*10+06)*coeffd(ilaux)
      xcfitd(3)=xcfit3(nxcs+(naux-1)*10+07)*coeffd(ilaux)
      xcfitd(4)=xcfit3(nxcs+(naux-1)*10+08)*coeffd(ilaux)/sqrt3
      xcfitd(5)=xcfit3(nxcs+(naux-1)*10+09)*coeffd(ilaux)
      xcfitd(6)=xcfit3(nxcs+(naux-1)*10+10)*coeffd(ilaux)/sqrt3

      do 2002 i=1,n
      zetafac(i)=(zeta(i)/(zeta(i)+zetaux))**(three/two)
 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 efac(i)=zetaux*zeta(i)*zetarec(i)*(pxcx*pxcx+pycy*pycy+pzcz*pzcz)

      do 2004 i=1,n
      sss(i)=zero
 2004 if(efac(i).lt.expcut) sss(i)=zetafac(i)*base(i)*exp(-efac(i))

      nn=0

      do 2005 i=1,n
      if(abs(sss(i)).gt.calcut) then
                                      nn=nn+1

                                      iloc(nn)=i

                                      iicontr(nn)=icontr(i)

                                      zzeta(nn)=zeta(i)

                                      zetarec(nn)=zetarec(i)

                                      aax(nn)=ax(i)
                                      aay(nn)=ay(i)
                                      aaz(nn)=az(i)

                                      axcx(nn)=ax(i)-ccx
                                      aycy(nn)=ay(i)-ccy
                                      azcz(nn)=az(i)-ccz

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

                                      sss(nn)=sss(i)
                                endif
 2005 continue

      do 2006 i=1,nn
      gxax(i)=(zzeta(i)*ppx(i)+zetaux*ccx)*zetarec(i)-aax(i)
      gyay(i)=(zzeta(i)*ppy(i)+zetaux*ccy)*zetarec(i)-aay(i)
 2006 gzaz(i)=(zzeta(i)*ppz(i)+zetaux*ccz)*zetarec(i)-aaz(i)

      do 2007 i=1,nn
 2007 zetarec(i)=zetarec(i)/two

      do 2008 i=1,nn
      pss(i,1)=gxax(i)*sss(i)
      pss(i,2)=gyay(i)*sss(i)
 2008 pss(i,3)=gzaz(i)*sss(i)

      do 2009 i=1,nn

      c(i)=zetarec(i)*sss(i)

      dss(i,1)=gxax(i)*pss(i,1)+c(i)
      dss(i,2)=gyay(i)*pss(i,1)
      dss(i,3)=gzaz(i)*pss(i,1)
      dss(i,4)=gyay(i)*pss(i,2)+c(i)
      dss(i,5)=gzaz(i)*pss(i,2)
      dss(i,6)=gzaz(i)*pss(i,3)+c(i)

 2009 continue

      do 2010 i=1,nn
      ssp(i,1)=axcx(i)*sss(i)+pss(i,1)
      ssp(i,2)=aycy(i)*sss(i)+pss(i,2)
 2010 ssp(i,3)=azcz(i)*sss(i)+pss(i,3)

      do 2011 i=1,nn
      psp(i,1)=axcx(i)*pss(i,1)+dss(i,1)
      psp(i,2)=aycy(i)*pss(i,1)+dss(i,2)
      psp(i,3)=azcz(i)*pss(i,1)+dss(i,3)
      psp(i,4)=axcx(i)*pss(i,2)+dss(i,2)
      psp(i,5)=aycy(i)*pss(i,2)+dss(i,4)
      psp(i,6)=azcz(i)*pss(i,2)+dss(i,5)
      psp(i,7)=axcx(i)*pss(i,3)+dss(i,3)
      psp(i,8)=aycy(i)*pss(i,3)+dss(i,5)
 2011 psp(i,9)=azcz(i)*pss(i,3)+dss(i,6)

      do 2012 i=1,nn
      ssd(i,1)=axcx(i)*ssp(i,1)+psp(i,1)
      ssd(i,2)=aycy(i)*ssp(i,1)+psp(i,4)
      ssd(i,3)=azcz(i)*ssp(i,1)+psp(i,7)
      ssd(i,4)=aycy(i)*ssp(i,2)+psp(i,5)
      ssd(i,5)=azcz(i)*ssp(i,2)+psp(i,8)
 2012 ssd(i,6)=azcz(i)*ssp(i,3)+psp(i,9)

cfpp$ nodepchk
      do 2013 i=1,nn
      index=iloc(i)
 2013 gather(index,1)=gather(index,1)+sss(i)*xcfits
     &                               +ssp(i,1)*xcfitp(1)
     &                               +ssp(i,2)*xcfitp(2)
     &                               +ssp(i,3)*xcfitp(3)
     &                               +ssd(i,1)*xcfitd(1)
     &                               +ssd(i,2)*xcfitd(2)
     &                               +ssd(i,3)*xcfitd(3)
     &                               +ssd(i,4)*xcfitd(4)
     &                               +ssd(i,5)*xcfitd(5)
     &                               +ssd(i,6)*xcfitd(6)

 2001 continue

      do 1006 i=1,n
 1006 fock3(inddyy(icontr(i)))=fock3(inddyy(icontr(i)))+gather(i,1)

      return
      end
