      subroutine xvecps(m,n,nxcs,nxcspd,ilfunc,icfunc,icontr,iicontr,
     &                  inddyy,iloc,coord,dmat,vvector,alphaaux,
     &                  coeffs,coeffp,coeffd,zeta,ax,ay,az,px,py,pz,
     &                  base,zetarec,zetafac,efac,gxax,gyay,gzaz,sss,
     &                  pss,psscond,dss,dsscond,fss,fsscond,psp,dsp,psd,
     &                  c,c1,zzeta,aax,aay,aaz,ppx,ppy,ppz,axcx,aycy,
     &                  azcz,gather)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
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(m,*),
     &          iloc(*)

      dimension coord(3,*),dmat(*),vvector(*),alphaaux(*),coeffs(*),
     &          coeffp(*),coeffd(*),zeta(*),ax(*),ay(*),az(*),px(*),
     &          py(*),pz(*),base(*),zetarec(*),zetafac(*),efac(*),
     &          gxax(*),gyay(*),gzaz(*),sss(*),pss(n,*),psscond(n,*),
     &          dss(n,*),dsscond(n,*),fss(n,*),fsscond(n,*),psp(n,*),
     &          dsp(n,*),psd(n,*),c(*),c1(n,*),zzeta(*),aax(*),aay(*),
     &          aaz(*),ppx(*),ppy(*),ppz(*),axcx(*),aycy(*),azcz(*),
     &          gather(n,*)

      data zero,one,two,three/0.0d0,1.0d0,2.0d0,3.0d0/
      data expcut,calcut,pi/45.0d0,1.0d-16,3.1415926535898d0/

      thrhlf=three/two

      sqrt3=sqrt(three)

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

      do 1001 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=coeffs(ilaux)

      do 1002 i=1,n
      zetafac(i)=(zeta(i)/(zeta(i)+zetaux))**(three/two)
 1002 zetarec(i)=one/(zeta(i)+zetaux)

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

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

      nn=0

      do 1005 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)

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

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

      do 1006 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)
 1006 gzaz(i)=(zzeta(i)*ppz(i)+zetaux*ccz)*zetarec(i)-aaz(i)

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

      do 1008 i=1,nn
 1008 gather(i,1)=zero

      do 1009 j=1,3
      do 1009 i=1,nn
 1009 gather(i,1)=gather(i,1)+dmat(inddyy(iicontr(i),j))*pss(i,j)

      do 1010 i=1,nn
 1010 gather(i,1)=gather(i,1)*xcfits

      do 1011 i=1,nn
 1011 vvector(naux)=vvector(naux)+gather(i,1)

 1001 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=coeffs(ilaux)

      xcfitp(1)=coeffp(ilaux)
      xcfitp(2)=coeffp(ilaux)
      xcfitp(3)=coeffp(ilaux)

      xcfitd(1)=coeffd(ilaux)/sqrt3
      xcfitd(2)=coeffd(ilaux)
      xcfitd(3)=coeffd(ilaux)
      xcfitd(4)=coeffd(ilaux)/sqrt3
      xcfitd(5)=coeffd(ilaux)
      xcfitd(6)=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 j=1,3
      do 2010 i=1,nn
 2010 c1(i,j)=zetarec(i)*pss(i,j)

      do 2011 i=1,nn
      fss(i,01)=gxax(i)*dss(i,1)+c1(i,1)+c1(i,1)
      fss(i,02)=gyay(i)*dss(i,1)
      fss(i,03)=gzaz(i)*dss(i,1)
      fss(i,04)=gyay(i)*dss(i,2)+c1(i,1)
      fss(i,05)=gzaz(i)*dss(i,2)
      fss(i,06)=gzaz(i)*dss(i,3)+c1(i,1)
      fss(i,07)=gyay(i)*dss(i,4)+c1(i,2)+c1(i,2)
      fss(i,08)=gzaz(i)*dss(i,4)
      fss(i,09)=gzaz(i)*dss(i,5)+c1(i,2)
 2011 fss(i,10)=gzaz(i)*dss(i,6)+c1(i,3)+c1(i,3)

      do 2012 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)
 2012 psp(i,9)=azcz(i)*pss(i,3)+dss(i,6)

      do 2013 i=1,nn
      dsp(i,01)=axcx(i)*dss(i,1)+fss(i,01)
      dsp(i,02)=aycy(i)*dss(i,1)+fss(i,02)
      dsp(i,03)=azcz(i)*dss(i,1)+fss(i,03)
      dsp(i,04)=axcx(i)*dss(i,2)+fss(i,02)
      dsp(i,05)=aycy(i)*dss(i,2)+fss(i,04)
      dsp(i,06)=azcz(i)*dss(i,2)+fss(i,05)
      dsp(i,07)=axcx(i)*dss(i,3)+fss(i,03)
      dsp(i,08)=aycy(i)*dss(i,3)+fss(i,05)
      dsp(i,09)=azcz(i)*dss(i,3)+fss(i,06)
      dsp(i,10)=axcx(i)*dss(i,4)+fss(i,04)
      dsp(i,11)=aycy(i)*dss(i,4)+fss(i,07)
      dsp(i,12)=azcz(i)*dss(i,4)+fss(i,08)
      dsp(i,13)=axcx(i)*dss(i,5)+fss(i,05)
      dsp(i,14)=aycy(i)*dss(i,5)+fss(i,08)
      dsp(i,15)=azcz(i)*dss(i,5)+fss(i,09)
      dsp(i,16)=axcx(i)*dss(i,6)+fss(i,06)
      dsp(i,17)=aycy(i)*dss(i,6)+fss(i,09)
 2013 dsp(i,18)=azcz(i)*dss(i,6)+fss(i,10)

      do 2014 i=1,nn
      psd(i,01)=axcx(i)*psp(i,1)+dsp(i,01)
      psd(i,02)=aycy(i)*psp(i,1)+dsp(i,04)
      psd(i,03)=azcz(i)*psp(i,1)+dsp(i,07)
      psd(i,04)=aycy(i)*psp(i,2)+dsp(i,05)
      psd(i,05)=azcz(i)*psp(i,2)+dsp(i,08)
      psd(i,06)=azcz(i)*psp(i,3)+dsp(i,09)
      psd(i,07)=axcx(i)*psp(i,4)+dsp(i,04)
      psd(i,08)=aycy(i)*psp(i,4)+dsp(i,10)
      psd(i,09)=azcz(i)*psp(i,4)+dsp(i,13)
      psd(i,10)=aycy(i)*psp(i,5)+dsp(i,11)
      psd(i,11)=azcz(i)*psp(i,5)+dsp(i,14)
      psd(i,12)=azcz(i)*psp(i,6)+dsp(i,15)
      psd(i,13)=axcx(i)*psp(i,7)+dsp(i,07)
      psd(i,14)=aycy(i)*psp(i,7)+dsp(i,13)
      psd(i,15)=azcz(i)*psp(i,7)+dsp(i,16)
      psd(i,16)=aycy(i)*psp(i,8)+dsp(i,14)
      psd(i,17)=azcz(i)*psp(i,8)+dsp(i,17)
 2014 psd(i,18)=azcz(i)*psp(i,9)+dsp(i,18)

      do 2015 j=1,10
      do 2015 i=1,nn
 2015 gather(i,j)=zero

      do 2016 j=1,3
      do 2016 i=1,nn
 2016 gather(i,1)=gather(i,1)+dmat(inddyy(iicontr(i),j))*pss(i,j)

      do 2017 k=1,3
      do 2017 j=1,3
      do 2017 i=1,nn
 2017 gather(i,k+1)=gather(i,k+1)
     &             +dmat(inddyy(iicontr(i),j))*psp(i,3*(j-1)+k)

      do 2018 k=1,6
      do 2018 j=1,3
      do 2018 i=1,nn
 2018 gather(i,k+4)=gather(i,k+4)
     &             +dmat(inddyy(iicontr(i),j))*psd(i,6*(j-1)+k)

      do 2019 i=1,nn
      gather(i,01)=gather(i,01)*xcfits
      gather(i,02)=gather(i,02)*xcfitp(1)
      gather(i,03)=gather(i,03)*xcfitp(2)
      gather(i,04)=gather(i,04)*xcfitp(3)
      gather(i,05)=gather(i,05)*xcfitd(1)
      gather(i,06)=gather(i,06)*xcfitd(2)
      gather(i,07)=gather(i,07)*xcfitd(3)
      gather(i,08)=gather(i,08)*xcfitd(4)
      gather(i,09)=gather(i,09)*xcfitd(5)
 2019 gather(i,10)=gather(i,10)*xcfitd(6)

      nc=nxcs+(naux-1)*10

      do 2020 i=1,nn
      do 2020 j=1,10
 2020 vvector(nc+j)=vvector(nc+j)+gather(i,j)

 2001 continue

      return
      end
