      subroutine xvecdp(m,n,nxcs,nxcspd,ilfunc,icfunc,icontr,iicontr,
     &                  inddyy,iloc,coord,dmat,vvector,alphaaux,
     &                  coeffs,coeffp,coeffd,zeta,ax,ay,az,axbx,ayby,
     &                  azbz,px,py,pz,base,zetarec,zetafac,efac,gxax,
     &                  gyay,gzaz,sss,pss,dss,dsscond,fss,fsscond,gss,
     &                  gsscond,hss,hsscond,dps,fps,gps,dpp,fpp,dpd,c,
     &                  c1,c2,c3,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(*),axbx(*),
     &          ayby(*),azbz(*),px(*),py(*),pz(*),base(*),zetarec(*),
     &          zetafac(*),efac(*),gxax(*),gyay(*),gzaz(*),sss(*),
     &          pss(n,*),dss(n,*),dsscond(n,*),fss(n,*),fsscond(n,*),
     &          gss(n,*),gsscond(n,*),hss(n,*),hsscond(n,*),dps(n,*),
     &          fps(n,*),gps(n,*),dpp(n,*),fpp(n,*),dpd(n,*),c(*),
     &          c1(n,*),c2(n,*),c3(n,*),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,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
 1007 zetarec(i)=zetarec(i)/two

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

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

 1009 continue

      do 1010 j=1,3
      do 1010 i=1,nn
 1010 c1(i,j)=zetarec(i)*pss(i,j)

      do 1011 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)
 1011 fss(i,10)=gzaz(i)*dss(i,6)+c1(i,3)+c1(i,3)

      do 1012 i=1,nn
      dps(i,01)=axbx(iicontr(i))*dss(i,1)+fss(i,01)
      dps(i,02)=ayby(iicontr(i))*dss(i,1)+fss(i,02)
      dps(i,03)=azbz(iicontr(i))*dss(i,1)+fss(i,03)
      dps(i,04)=axbx(iicontr(i))*dss(i,2)+fss(i,02)
      dps(i,05)=ayby(iicontr(i))*dss(i,2)+fss(i,04)
      dps(i,06)=azbz(iicontr(i))*dss(i,2)+fss(i,05)
      dps(i,07)=axbx(iicontr(i))*dss(i,3)+fss(i,03)
      dps(i,08)=ayby(iicontr(i))*dss(i,3)+fss(i,05)
      dps(i,09)=azbz(iicontr(i))*dss(i,3)+fss(i,06)
      dps(i,10)=axbx(iicontr(i))*dss(i,4)+fss(i,04)
      dps(i,11)=ayby(iicontr(i))*dss(i,4)+fss(i,07)
      dps(i,12)=azbz(iicontr(i))*dss(i,4)+fss(i,08)
      dps(i,13)=axbx(iicontr(i))*dss(i,5)+fss(i,05)
      dps(i,14)=ayby(iicontr(i))*dss(i,5)+fss(i,08)
      dps(i,15)=azbz(iicontr(i))*dss(i,5)+fss(i,09)
      dps(i,16)=axbx(iicontr(i))*dss(i,6)+fss(i,06)
      dps(i,17)=ayby(iicontr(i))*dss(i,6)+fss(i,09)
 1012 dps(i,18)=azbz(iicontr(i))*dss(i,6)+fss(i,10)

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

      do 1014 j=1,18
      do 1014 i=1,nn
 1014 gather(i,1)=gather(i,1)+dmat(inddyy(iicontr(i),j))*dps(i,j)

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

      do 1016 i=1,nn
 1016 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 j=1,6
      do 2012 i=1,nn
 2012 c2(i,j)=zetarec(i)*dss(i,j)

      do 2013 i=1,nn
      gss(i,01)=gxax(i)*fss(i,01)+c2(i,1)*three
      gss(i,02)=gyay(i)*fss(i,01)
      gss(i,03)=gzaz(i)*fss(i,01)
      gss(i,04)=gyay(i)*fss(i,02)+c2(i,1)
      gss(i,05)=gzaz(i)*fss(i,02)
      gss(i,06)=gzaz(i)*fss(i,03)+c2(i,1)
      gss(i,07)=gyay(i)*fss(i,04)+c2(i,2)*two
      gss(i,08)=gzaz(i)*fss(i,04)
      gss(i,09)=gzaz(i)*fss(i,05)+c2(i,2)
      gss(i,10)=gzaz(i)*fss(i,06)+c2(i,3)*two
      gss(i,11)=gyay(i)*fss(i,07)+c2(i,4)*three
      gss(i,12)=gzaz(i)*fss(i,07)
      gss(i,13)=gzaz(i)*fss(i,08)+c2(i,4)
      gss(i,14)=gzaz(i)*fss(i,09)+c2(i,5)*two
 2013 gss(i,15)=gzaz(i)*fss(i,10)+c2(i,6)*three

      do 2014 j=1,10
      do 2014 i=1,nn
 2014 c3(i,j)=zetarec(i)*fss(i,j)

      do 2015 i=1,nn
      hss(i,01)=gxax(i)*gss(i,01)+c3(i,01)*four
      hss(i,02)=gyay(i)*gss(i,01)
      hss(i,03)=gzaz(i)*gss(i,01)
      hss(i,04)=gyay(i)*gss(i,02)+c3(i,01)
      hss(i,05)=gzaz(i)*gss(i,02)
      hss(i,06)=gzaz(i)*gss(i,03)+c3(i,01)
      hss(i,07)=gyay(i)*gss(i,04)+c3(i,02)*two
      hss(i,08)=gzaz(i)*gss(i,04)
      hss(i,09)=gzaz(i)*gss(i,05)+c3(i,02)
      hss(i,10)=gzaz(i)*gss(i,06)+c3(i,03)*two
      hss(i,11)=gyay(i)*gss(i,07)+c3(i,04)*three
      hss(i,12)=gzaz(i)*gss(i,07)
      hss(i,13)=gzaz(i)*gss(i,08)+c3(i,04)
      hss(i,14)=gzaz(i)*gss(i,09)+c3(i,05)*two
      hss(i,15)=gzaz(i)*gss(i,10)+c3(i,06)*three
      hss(i,16)=gyay(i)*gss(i,11)+c3(i,07)*four
      hss(i,17)=gzaz(i)*gss(i,11)
      hss(i,18)=gzaz(i)*gss(i,12)+c3(i,07)
      hss(i,19)=gzaz(i)*gss(i,13)+c3(i,08)*two
      hss(i,20)=gzaz(i)*gss(i,14)+c3(i,09)*three
 2015 hss(i,21)=gzaz(i)*gss(i,15)+c3(i,10)*four

      do 2016 i=1,nn
      dps(i,01)=axbx(iicontr(i))*dss(i,1)+fss(i,01)
      dps(i,02)=ayby(iicontr(i))*dss(i,1)+fss(i,02)
      dps(i,03)=azbz(iicontr(i))*dss(i,1)+fss(i,03)
      dps(i,04)=axbx(iicontr(i))*dss(i,2)+fss(i,02)
      dps(i,05)=ayby(iicontr(i))*dss(i,2)+fss(i,04)
      dps(i,06)=azbz(iicontr(i))*dss(i,2)+fss(i,05)
      dps(i,07)=axbx(iicontr(i))*dss(i,3)+fss(i,03)
      dps(i,08)=ayby(iicontr(i))*dss(i,3)+fss(i,05)
      dps(i,09)=azbz(iicontr(i))*dss(i,3)+fss(i,06)
      dps(i,10)=axbx(iicontr(i))*dss(i,4)+fss(i,04)
      dps(i,11)=ayby(iicontr(i))*dss(i,4)+fss(i,07)
      dps(i,12)=azbz(iicontr(i))*dss(i,4)+fss(i,08)
      dps(i,13)=axbx(iicontr(i))*dss(i,5)+fss(i,05)
      dps(i,14)=ayby(iicontr(i))*dss(i,5)+fss(i,08)
      dps(i,15)=azbz(iicontr(i))*dss(i,5)+fss(i,09)
      dps(i,16)=axbx(iicontr(i))*dss(i,6)+fss(i,06)
      dps(i,17)=ayby(iicontr(i))*dss(i,6)+fss(i,09)
 2016 dps(i,18)=azbz(iicontr(i))*dss(i,6)+fss(i,10)

      do 2017 i=1,nn
      fps(i,01)=axbx(iicontr(i))*fss(i,01)+gss(i,01)
      fps(i,02)=ayby(iicontr(i))*fss(i,01)+gss(i,02)
      fps(i,03)=azbz(iicontr(i))*fss(i,01)+gss(i,03)
      fps(i,04)=axbx(iicontr(i))*fss(i,02)+gss(i,02)
      fps(i,05)=ayby(iicontr(i))*fss(i,02)+gss(i,04)
      fps(i,06)=azbz(iicontr(i))*fss(i,02)+gss(i,05)
      fps(i,07)=axbx(iicontr(i))*fss(i,03)+gss(i,03)
      fps(i,08)=ayby(iicontr(i))*fss(i,03)+gss(i,05)
      fps(i,09)=azbz(iicontr(i))*fss(i,03)+gss(i,06)
      fps(i,10)=axbx(iicontr(i))*fss(i,04)+gss(i,04)
      fps(i,11)=ayby(iicontr(i))*fss(i,04)+gss(i,07)
      fps(i,12)=azbz(iicontr(i))*fss(i,04)+gss(i,08)
      fps(i,13)=axbx(iicontr(i))*fss(i,05)+gss(i,05)
      fps(i,14)=ayby(iicontr(i))*fss(i,05)+gss(i,08)
      fps(i,15)=azbz(iicontr(i))*fss(i,05)+gss(i,09)
      fps(i,16)=axbx(iicontr(i))*fss(i,06)+gss(i,06)
      fps(i,17)=ayby(iicontr(i))*fss(i,06)+gss(i,09)
      fps(i,18)=azbz(iicontr(i))*fss(i,06)+gss(i,10)
      fps(i,19)=axbx(iicontr(i))*fss(i,07)+gss(i,07)
      fps(i,20)=ayby(iicontr(i))*fss(i,07)+gss(i,11)
      fps(i,21)=azbz(iicontr(i))*fss(i,07)+gss(i,12)
      fps(i,22)=axbx(iicontr(i))*fss(i,08)+gss(i,08)
      fps(i,23)=ayby(iicontr(i))*fss(i,08)+gss(i,12)
      fps(i,24)=azbz(iicontr(i))*fss(i,08)+gss(i,13)
      fps(i,25)=axbx(iicontr(i))*fss(i,09)+gss(i,09)
      fps(i,26)=ayby(iicontr(i))*fss(i,09)+gss(i,13)
      fps(i,27)=azbz(iicontr(i))*fss(i,09)+gss(i,14)
      fps(i,28)=axbx(iicontr(i))*fss(i,10)+gss(i,10)
      fps(i,29)=ayby(iicontr(i))*fss(i,10)+gss(i,14)
 2017 fps(i,30)=azbz(iicontr(i))*fss(i,10)+gss(i,15)

      do 2018 i=1,nn
      dpp(i,01)=axcx(i)*dps(i,01)+fps(i,01)
      dpp(i,02)=aycy(i)*dps(i,01)+fps(i,04)
      dpp(i,03)=azcz(i)*dps(i,01)+fps(i,07)
      dpp(i,04)=axcx(i)*dps(i,02)+fps(i,02)
      dpp(i,05)=aycy(i)*dps(i,02)+fps(i,05)
      dpp(i,06)=azcz(i)*dps(i,02)+fps(i,08)
      dpp(i,07)=axcx(i)*dps(i,03)+fps(i,03)
      dpp(i,08)=aycy(i)*dps(i,03)+fps(i,06)
      dpp(i,09)=azcz(i)*dps(i,03)+fps(i,09)
      dpp(i,10)=axcx(i)*dps(i,04)+fps(i,04)
      dpp(i,11)=aycy(i)*dps(i,04)+fps(i,10)
      dpp(i,12)=azcz(i)*dps(i,04)+fps(i,13)
      dpp(i,13)=axcx(i)*dps(i,05)+fps(i,05)
      dpp(i,14)=aycy(i)*dps(i,05)+fps(i,11)
      dpp(i,15)=azcz(i)*dps(i,05)+fps(i,14)
      dpp(i,16)=axcx(i)*dps(i,06)+fps(i,06)
      dpp(i,17)=aycy(i)*dps(i,06)+fps(i,12)
      dpp(i,18)=azcz(i)*dps(i,06)+fps(i,15)
      dpp(i,19)=axcx(i)*dps(i,07)+fps(i,07)
      dpp(i,20)=aycy(i)*dps(i,07)+fps(i,13)
      dpp(i,21)=azcz(i)*dps(i,07)+fps(i,16)
      dpp(i,22)=axcx(i)*dps(i,08)+fps(i,08)
      dpp(i,23)=aycy(i)*dps(i,08)+fps(i,14)
      dpp(i,24)=azcz(i)*dps(i,08)+fps(i,17)
      dpp(i,25)=axcx(i)*dps(i,09)+fps(i,09)
      dpp(i,26)=aycy(i)*dps(i,09)+fps(i,15)
      dpp(i,27)=azcz(i)*dps(i,09)+fps(i,18)
      dpp(i,28)=axcx(i)*dps(i,10)+fps(i,10)
      dpp(i,29)=aycy(i)*dps(i,10)+fps(i,19)
      dpp(i,30)=azcz(i)*dps(i,10)+fps(i,22)
      dpp(i,31)=axcx(i)*dps(i,11)+fps(i,11)
      dpp(i,32)=aycy(i)*dps(i,11)+fps(i,20)
      dpp(i,33)=azcz(i)*dps(i,11)+fps(i,23)
      dpp(i,34)=axcx(i)*dps(i,12)+fps(i,12)
      dpp(i,35)=aycy(i)*dps(i,12)+fps(i,21)
      dpp(i,36)=azcz(i)*dps(i,12)+fps(i,24)
      dpp(i,37)=axcx(i)*dps(i,13)+fps(i,13)
      dpp(i,38)=aycy(i)*dps(i,13)+fps(i,22)
      dpp(i,39)=azcz(i)*dps(i,13)+fps(i,25)
      dpp(i,40)=axcx(i)*dps(i,14)+fps(i,14)
      dpp(i,41)=aycy(i)*dps(i,14)+fps(i,23)
      dpp(i,42)=azcz(i)*dps(i,14)+fps(i,26)
      dpp(i,43)=axcx(i)*dps(i,15)+fps(i,15)
      dpp(i,44)=aycy(i)*dps(i,15)+fps(i,24)
      dpp(i,45)=azcz(i)*dps(i,15)+fps(i,27)
      dpp(i,46)=axcx(i)*dps(i,16)+fps(i,16)
      dpp(i,47)=aycy(i)*dps(i,16)+fps(i,25)
      dpp(i,48)=azcz(i)*dps(i,16)+fps(i,28)
      dpp(i,49)=axcx(i)*dps(i,17)+fps(i,17)
      dpp(i,50)=aycy(i)*dps(i,17)+fps(i,26)
      dpp(i,51)=azcz(i)*dps(i,17)+fps(i,29)
      dpp(i,52)=axcx(i)*dps(i,18)+fps(i,18)
      dpp(i,53)=aycy(i)*dps(i,18)+fps(i,27)
 2018 dpp(i,54)=azcz(i)*dps(i,18)+fps(i,30)

      do 2019 i=1,nn
      gps(i,01)=axbx(iicontr(i))*gss(i,01)+hss(i,01)
      gps(i,02)=ayby(iicontr(i))*gss(i,01)+hss(i,02)
      gps(i,03)=azbz(iicontr(i))*gss(i,01)+hss(i,03)
      gps(i,04)=axbx(iicontr(i))*gss(i,02)+hss(i,02)
      gps(i,05)=ayby(iicontr(i))*gss(i,02)+hss(i,04)
      gps(i,06)=azbz(iicontr(i))*gss(i,02)+hss(i,05)
      gps(i,07)=axbx(iicontr(i))*gss(i,03)+hss(i,03)
      gps(i,08)=ayby(iicontr(i))*gss(i,03)+hss(i,05)
      gps(i,09)=azbz(iicontr(i))*gss(i,03)+hss(i,06)
      gps(i,10)=axbx(iicontr(i))*gss(i,04)+hss(i,04)
      gps(i,11)=ayby(iicontr(i))*gss(i,04)+hss(i,07)
      gps(i,12)=azbz(iicontr(i))*gss(i,04)+hss(i,08)
      gps(i,13)=axbx(iicontr(i))*gss(i,05)+hss(i,05)
      gps(i,14)=ayby(iicontr(i))*gss(i,05)+hss(i,08)
      gps(i,15)=azbz(iicontr(i))*gss(i,05)+hss(i,09)
      gps(i,16)=axbx(iicontr(i))*gss(i,06)+hss(i,06)
      gps(i,17)=ayby(iicontr(i))*gss(i,06)+hss(i,09)
      gps(i,18)=azbz(iicontr(i))*gss(i,06)+hss(i,10)
      gps(i,19)=axbx(iicontr(i))*gss(i,07)+hss(i,07)
      gps(i,20)=ayby(iicontr(i))*gss(i,07)+hss(i,11)
      gps(i,21)=azbz(iicontr(i))*gss(i,07)+hss(i,12)
      gps(i,22)=axbx(iicontr(i))*gss(i,08)+hss(i,08)
      gps(i,23)=ayby(iicontr(i))*gss(i,08)+hss(i,12)
      gps(i,24)=azbz(iicontr(i))*gss(i,08)+hss(i,13)
      gps(i,25)=axbx(iicontr(i))*gss(i,09)+hss(i,09)
      gps(i,26)=ayby(iicontr(i))*gss(i,09)+hss(i,13)
      gps(i,27)=azbz(iicontr(i))*gss(i,09)+hss(i,14)
      gps(i,28)=axbx(iicontr(i))*gss(i,10)+hss(i,10)
      gps(i,29)=ayby(iicontr(i))*gss(i,10)+hss(i,14)
      gps(i,30)=azbz(iicontr(i))*gss(i,10)+hss(i,15)
      gps(i,31)=axbx(iicontr(i))*gss(i,11)+hss(i,11)
      gps(i,32)=ayby(iicontr(i))*gss(i,11)+hss(i,16)
      gps(i,33)=azbz(iicontr(i))*gss(i,11)+hss(i,17)
      gps(i,34)=axbx(iicontr(i))*gss(i,12)+hss(i,12)
      gps(i,35)=ayby(iicontr(i))*gss(i,12)+hss(i,17)
      gps(i,36)=azbz(iicontr(i))*gss(i,12)+hss(i,18)
      gps(i,37)=axbx(iicontr(i))*gss(i,13)+hss(i,13)
      gps(i,38)=ayby(iicontr(i))*gss(i,13)+hss(i,18)
      gps(i,39)=azbz(iicontr(i))*gss(i,13)+hss(i,19)
      gps(i,40)=axbx(iicontr(i))*gss(i,14)+hss(i,14)
      gps(i,41)=ayby(iicontr(i))*gss(i,14)+hss(i,19)
      gps(i,42)=azbz(iicontr(i))*gss(i,14)+hss(i,20)
      gps(i,43)=axbx(iicontr(i))*gss(i,15)+hss(i,15)
      gps(i,44)=ayby(iicontr(i))*gss(i,15)+hss(i,20)
 2019 gps(i,45)=azbz(iicontr(i))*gss(i,15)+hss(i,21)

      do 2020 i=1,nn
      fpp(i,01)=axcx(i)*fps(i,01)+gps(i,01)
      fpp(i,02)=aycy(i)*fps(i,01)+gps(i,04)
      fpp(i,03)=azcz(i)*fps(i,01)+gps(i,07)
      fpp(i,04)=axcx(i)*fps(i,02)+gps(i,02)
      fpp(i,05)=aycy(i)*fps(i,02)+gps(i,05)
      fpp(i,06)=azcz(i)*fps(i,02)+gps(i,08)
      fpp(i,07)=axcx(i)*fps(i,03)+gps(i,03)
      fpp(i,08)=aycy(i)*fps(i,03)+gps(i,06)
      fpp(i,09)=azcz(i)*fps(i,03)+gps(i,09)
      fpp(i,10)=axcx(i)*fps(i,04)+gps(i,04)
      fpp(i,11)=aycy(i)*fps(i,04)+gps(i,10)
      fpp(i,12)=azcz(i)*fps(i,04)+gps(i,13)
      fpp(i,13)=axcx(i)*fps(i,05)+gps(i,05)
      fpp(i,14)=aycy(i)*fps(i,05)+gps(i,11)
      fpp(i,15)=azcz(i)*fps(i,05)+gps(i,14)
      fpp(i,16)=axcx(i)*fps(i,06)+gps(i,06)
      fpp(i,17)=aycy(i)*fps(i,06)+gps(i,12)
      fpp(i,18)=azcz(i)*fps(i,06)+gps(i,15)
      fpp(i,19)=axcx(i)*fps(i,07)+gps(i,07)
      fpp(i,20)=aycy(i)*fps(i,07)+gps(i,13)
      fpp(i,21)=azcz(i)*fps(i,07)+gps(i,16)
      fpp(i,22)=axcx(i)*fps(i,08)+gps(i,08)
      fpp(i,23)=aycy(i)*fps(i,08)+gps(i,14)
      fpp(i,24)=azcz(i)*fps(i,08)+gps(i,17)
      fpp(i,25)=axcx(i)*fps(i,09)+gps(i,09)
      fpp(i,26)=aycy(i)*fps(i,09)+gps(i,15)
      fpp(i,27)=azcz(i)*fps(i,09)+gps(i,18)
      fpp(i,28)=axcx(i)*fps(i,10)+gps(i,10)
      fpp(i,29)=aycy(i)*fps(i,10)+gps(i,19)
      fpp(i,30)=azcz(i)*fps(i,10)+gps(i,22)
      fpp(i,31)=axcx(i)*fps(i,11)+gps(i,11)
      fpp(i,32)=aycy(i)*fps(i,11)+gps(i,20)
      fpp(i,33)=azcz(i)*fps(i,11)+gps(i,23)
      fpp(i,34)=axcx(i)*fps(i,12)+gps(i,12)
      fpp(i,35)=aycy(i)*fps(i,12)+gps(i,21)
      fpp(i,36)=azcz(i)*fps(i,12)+gps(i,24)
      fpp(i,37)=axcx(i)*fps(i,13)+gps(i,13)
      fpp(i,38)=aycy(i)*fps(i,13)+gps(i,22)
      fpp(i,39)=azcz(i)*fps(i,13)+gps(i,25)
      fpp(i,40)=axcx(i)*fps(i,14)+gps(i,14)
      fpp(i,41)=aycy(i)*fps(i,14)+gps(i,23)
      fpp(i,42)=azcz(i)*fps(i,14)+gps(i,26)
      fpp(i,43)=axcx(i)*fps(i,15)+gps(i,15)
      fpp(i,44)=aycy(i)*fps(i,15)+gps(i,24)
      fpp(i,45)=azcz(i)*fps(i,15)+gps(i,27)
      fpp(i,46)=axcx(i)*fps(i,16)+gps(i,16)
      fpp(i,47)=aycy(i)*fps(i,16)+gps(i,25)
      fpp(i,48)=azcz(i)*fps(i,16)+gps(i,28)
      fpp(i,49)=axcx(i)*fps(i,17)+gps(i,17)
      fpp(i,50)=aycy(i)*fps(i,17)+gps(i,26)
      fpp(i,51)=azcz(i)*fps(i,17)+gps(i,29)
      fpp(i,52)=axcx(i)*fps(i,18)+gps(i,18)
      fpp(i,53)=aycy(i)*fps(i,18)+gps(i,27)
      fpp(i,54)=azcz(i)*fps(i,18)+gps(i,30)
      fpp(i,55)=axcx(i)*fps(i,19)+gps(i,19)
      fpp(i,56)=aycy(i)*fps(i,19)+gps(i,31)
      fpp(i,57)=azcz(i)*fps(i,19)+gps(i,34)
      fpp(i,58)=axcx(i)*fps(i,20)+gps(i,20)
      fpp(i,59)=aycy(i)*fps(i,20)+gps(i,32)
      fpp(i,60)=azcz(i)*fps(i,20)+gps(i,35)
      fpp(i,61)=axcx(i)*fps(i,21)+gps(i,21)
      fpp(i,62)=aycy(i)*fps(i,21)+gps(i,33)
      fpp(i,63)=azcz(i)*fps(i,21)+gps(i,36)
      fpp(i,64)=axcx(i)*fps(i,22)+gps(i,22)
      fpp(i,65)=aycy(i)*fps(i,22)+gps(i,34)
      fpp(i,66)=azcz(i)*fps(i,22)+gps(i,37)
      fpp(i,67)=axcx(i)*fps(i,23)+gps(i,23)
      fpp(i,68)=aycy(i)*fps(i,23)+gps(i,35)
      fpp(i,69)=azcz(i)*fps(i,23)+gps(i,38)
      fpp(i,70)=axcx(i)*fps(i,24)+gps(i,24)
      fpp(i,71)=aycy(i)*fps(i,24)+gps(i,36)
      fpp(i,72)=azcz(i)*fps(i,24)+gps(i,39)
      fpp(i,73)=axcx(i)*fps(i,25)+gps(i,25)
      fpp(i,74)=aycy(i)*fps(i,25)+gps(i,37)
      fpp(i,75)=azcz(i)*fps(i,25)+gps(i,40)
      fpp(i,76)=axcx(i)*fps(i,26)+gps(i,26)
      fpp(i,77)=aycy(i)*fps(i,26)+gps(i,38)
      fpp(i,78)=azcz(i)*fps(i,26)+gps(i,41)
      fpp(i,79)=axcx(i)*fps(i,27)+gps(i,27)
      fpp(i,80)=aycy(i)*fps(i,27)+gps(i,39)
      fpp(i,81)=azcz(i)*fps(i,27)+gps(i,42)
      fpp(i,82)=axcx(i)*fps(i,28)+gps(i,28)
      fpp(i,83)=aycy(i)*fps(i,28)+gps(i,40)
      fpp(i,84)=azcz(i)*fps(i,28)+gps(i,43)
      fpp(i,85)=axcx(i)*fps(i,29)+gps(i,29)
      fpp(i,86)=aycy(i)*fps(i,29)+gps(i,41)
      fpp(i,87)=azcz(i)*fps(i,29)+gps(i,44)
      fpp(i,88)=axcx(i)*fps(i,30)+gps(i,30)
      fpp(i,89)=aycy(i)*fps(i,30)+gps(i,42)
 2020 fpp(i,90)=azcz(i)*fps(i,30)+gps(i,45)

      do 2021 i=1,nn
      dpd(i,001)=axcx(i)*dpp(i,01)+fpp(i,01)
      dpd(i,002)=aycy(i)*dpp(i,01)+fpp(i,10)
      dpd(i,003)=azcz(i)*dpp(i,01)+fpp(i,19)
      dpd(i,004)=aycy(i)*dpp(i,02)+fpp(i,11)
      dpd(i,005)=azcz(i)*dpp(i,02)+fpp(i,20)
      dpd(i,006)=azcz(i)*dpp(i,03)+fpp(i,21)
      dpd(i,007)=axcx(i)*dpp(i,04)+fpp(i,04)
      dpd(i,008)=aycy(i)*dpp(i,04)+fpp(i,13)
      dpd(i,009)=azcz(i)*dpp(i,04)+fpp(i,22)
      dpd(i,010)=aycy(i)*dpp(i,05)+fpp(i,14)
      dpd(i,011)=azcz(i)*dpp(i,05)+fpp(i,23)
      dpd(i,012)=azcz(i)*dpp(i,06)+fpp(i,24)
      dpd(i,013)=axcx(i)*dpp(i,07)+fpp(i,07)
      dpd(i,014)=aycy(i)*dpp(i,07)+fpp(i,16)
      dpd(i,015)=azcz(i)*dpp(i,07)+fpp(i,25)
      dpd(i,016)=aycy(i)*dpp(i,08)+fpp(i,17)
      dpd(i,017)=azcz(i)*dpp(i,08)+fpp(i,26)
      dpd(i,018)=azcz(i)*dpp(i,09)+fpp(i,27)
      dpd(i,019)=axcx(i)*dpp(i,10)+fpp(i,10)
      dpd(i,020)=aycy(i)*dpp(i,10)+fpp(i,28)
      dpd(i,021)=azcz(i)*dpp(i,10)+fpp(i,37)
      dpd(i,022)=aycy(i)*dpp(i,11)+fpp(i,29)
      dpd(i,023)=azcz(i)*dpp(i,11)+fpp(i,38)
      dpd(i,024)=azcz(i)*dpp(i,12)+fpp(i,39)
      dpd(i,025)=axcx(i)*dpp(i,13)+fpp(i,13)
      dpd(i,026)=aycy(i)*dpp(i,13)+fpp(i,31)
      dpd(i,027)=azcz(i)*dpp(i,13)+fpp(i,40)
      dpd(i,028)=aycy(i)*dpp(i,14)+fpp(i,32)
      dpd(i,029)=azcz(i)*dpp(i,14)+fpp(i,41)
      dpd(i,030)=azcz(i)*dpp(i,15)+fpp(i,42)
      dpd(i,031)=axcx(i)*dpp(i,16)+fpp(i,16)
      dpd(i,032)=aycy(i)*dpp(i,16)+fpp(i,34)
      dpd(i,033)=azcz(i)*dpp(i,16)+fpp(i,43)
      dpd(i,034)=aycy(i)*dpp(i,17)+fpp(i,35)
      dpd(i,035)=azcz(i)*dpp(i,17)+fpp(i,44)
      dpd(i,036)=azcz(i)*dpp(i,18)+fpp(i,45)
      dpd(i,037)=axcx(i)*dpp(i,19)+fpp(i,19)
      dpd(i,038)=aycy(i)*dpp(i,19)+fpp(i,37)
      dpd(i,039)=azcz(i)*dpp(i,19)+fpp(i,46)
      dpd(i,040)=aycy(i)*dpp(i,20)+fpp(i,38)
      dpd(i,041)=azcz(i)*dpp(i,20)+fpp(i,47)
      dpd(i,042)=azcz(i)*dpp(i,21)+fpp(i,48)
      dpd(i,043)=axcx(i)*dpp(i,22)+fpp(i,22)
      dpd(i,044)=aycy(i)*dpp(i,22)+fpp(i,40)
      dpd(i,045)=azcz(i)*dpp(i,22)+fpp(i,49)
      dpd(i,046)=aycy(i)*dpp(i,23)+fpp(i,41)
      dpd(i,047)=azcz(i)*dpp(i,23)+fpp(i,50)
      dpd(i,048)=azcz(i)*dpp(i,24)+fpp(i,51)
      dpd(i,049)=axcx(i)*dpp(i,25)+fpp(i,25)
      dpd(i,050)=aycy(i)*dpp(i,25)+fpp(i,43)
      dpd(i,051)=azcz(i)*dpp(i,25)+fpp(i,52)
      dpd(i,052)=aycy(i)*dpp(i,26)+fpp(i,44)
      dpd(i,053)=azcz(i)*dpp(i,26)+fpp(i,53)
      dpd(i,054)=azcz(i)*dpp(i,27)+fpp(i,54)
      dpd(i,055)=axcx(i)*dpp(i,28)+fpp(i,28)
      dpd(i,056)=aycy(i)*dpp(i,28)+fpp(i,55)
      dpd(i,057)=azcz(i)*dpp(i,28)+fpp(i,64)
      dpd(i,058)=aycy(i)*dpp(i,29)+fpp(i,56)
      dpd(i,059)=azcz(i)*dpp(i,29)+fpp(i,65)
      dpd(i,060)=azcz(i)*dpp(i,30)+fpp(i,66)
      dpd(i,061)=axcx(i)*dpp(i,31)+fpp(i,31)
      dpd(i,062)=aycy(i)*dpp(i,31)+fpp(i,58)
      dpd(i,063)=azcz(i)*dpp(i,31)+fpp(i,67)
      dpd(i,064)=aycy(i)*dpp(i,32)+fpp(i,59)
      dpd(i,065)=azcz(i)*dpp(i,32)+fpp(i,68)
      dpd(i,066)=azcz(i)*dpp(i,33)+fpp(i,69)
      dpd(i,067)=axcx(i)*dpp(i,34)+fpp(i,34)
      dpd(i,068)=aycy(i)*dpp(i,34)+fpp(i,61)
      dpd(i,069)=azcz(i)*dpp(i,34)+fpp(i,70)
      dpd(i,070)=aycy(i)*dpp(i,35)+fpp(i,62)
      dpd(i,071)=azcz(i)*dpp(i,35)+fpp(i,71)
      dpd(i,072)=azcz(i)*dpp(i,36)+fpp(i,72)
      dpd(i,073)=axcx(i)*dpp(i,37)+fpp(i,37)
      dpd(i,074)=aycy(i)*dpp(i,37)+fpp(i,64)
      dpd(i,075)=azcz(i)*dpp(i,37)+fpp(i,73)
      dpd(i,076)=aycy(i)*dpp(i,38)+fpp(i,65)
      dpd(i,077)=azcz(i)*dpp(i,38)+fpp(i,74)
      dpd(i,078)=azcz(i)*dpp(i,39)+fpp(i,75)
      dpd(i,079)=axcx(i)*dpp(i,40)+fpp(i,40)
      dpd(i,080)=aycy(i)*dpp(i,40)+fpp(i,67)
      dpd(i,081)=azcz(i)*dpp(i,40)+fpp(i,76)
      dpd(i,082)=aycy(i)*dpp(i,41)+fpp(i,68)
      dpd(i,083)=azcz(i)*dpp(i,41)+fpp(i,77)
      dpd(i,084)=azcz(i)*dpp(i,42)+fpp(i,78)
      dpd(i,085)=axcx(i)*dpp(i,43)+fpp(i,43)
      dpd(i,086)=aycy(i)*dpp(i,43)+fpp(i,70)
      dpd(i,087)=azcz(i)*dpp(i,43)+fpp(i,79)
      dpd(i,088)=aycy(i)*dpp(i,44)+fpp(i,71)
      dpd(i,089)=azcz(i)*dpp(i,44)+fpp(i,80)
      dpd(i,090)=azcz(i)*dpp(i,45)+fpp(i,81)
      dpd(i,091)=axcx(i)*dpp(i,46)+fpp(i,46)
      dpd(i,092)=aycy(i)*dpp(i,46)+fpp(i,73)
      dpd(i,093)=azcz(i)*dpp(i,46)+fpp(i,82)
      dpd(i,094)=aycy(i)*dpp(i,47)+fpp(i,74)
      dpd(i,095)=azcz(i)*dpp(i,47)+fpp(i,83)
      dpd(i,096)=azcz(i)*dpp(i,48)+fpp(i,84)
      dpd(i,097)=axcx(i)*dpp(i,49)+fpp(i,49)
      dpd(i,098)=aycy(i)*dpp(i,49)+fpp(i,76)
      dpd(i,099)=azcz(i)*dpp(i,49)+fpp(i,85)
      dpd(i,100)=aycy(i)*dpp(i,50)+fpp(i,77)
      dpd(i,101)=azcz(i)*dpp(i,50)+fpp(i,86)
      dpd(i,102)=azcz(i)*dpp(i,51)+fpp(i,87)
      dpd(i,103)=axcx(i)*dpp(i,52)+fpp(i,52)
      dpd(i,104)=aycy(i)*dpp(i,52)+fpp(i,79)
      dpd(i,105)=azcz(i)*dpp(i,52)+fpp(i,88)
      dpd(i,106)=aycy(i)*dpp(i,53)+fpp(i,80)
      dpd(i,107)=azcz(i)*dpp(i,53)+fpp(i,89)
 2021 dpd(i,108)=azcz(i)*dpp(i,54)+fpp(i,90)

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

      do 2023 j=1,18
      do 2023 i=1,nn
 2023 gather(i,1)=gather(i,1)+dmat(inddyy(iicontr(i),j))*dps(i,j)

      do 2024 k=1,3
      do 2024 j=1,18
      do 2024 i=1,nn
 2024 gather(i,k+1)=gather(i,k+1)
     &             +dmat(inddyy(iicontr(i),j))*dpp(i,3*(j-1)+k)

      do 2025 k=1,6
      do 2025 j=1,18
      do 2025 i=1,nn
 2025 gather(i,k+4)=gather(i,k+4)
     &             +dmat(inddyy(iicontr(i),j))*dpd(i,6*(j-1)+k)

      do 2026 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)
 2026 gather(i,10)=gather(i,10)*xcfitd(6)

      nc=nxcs+(naux-1)*10

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

 2001 continue

      return
      end
