      subroutine xcordd(m,n,nxcs,nxcspd,ilfunc,icfunc,iminj,icontr,
     &                  iicontr,inddyy,iloc,coord,xcfit1,xcfit2,xcfit3,
     &                  alphaaux,coeffs,coeffp,coeffd,fock1,fock2,fock3,
     &                  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,qss,
     &                  qsscond,dps,dds,ddp,ddd,fps,fds,fdp,gps,gds,hps,
     &                  c,c1,c2,c3,c4,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 index(36)

      dimension xcfitp(3),xcfitd(6)

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

      dimension coord(3,*),xcfit1(*),xcfit2(*),xcfit3(*),alphaaux(*),
     &          coeffs(*),coeffp(*),coeffd(*),fock1(*),fock2(*),
     &          fock3(*),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,*),qss(n,*),
     &          qsscond(n,*),dps(n,*),dds(n,*),ddp(n,*),ddd(n,*),
     &          fps(n,*),fds(n,*),fdp(n,*),gps(n,*),gds(n,*),hps(n,*),
     &          c(*),c1(n,*),c2(n,*),c3(n,*),c4(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/45.0d0,1.0d-16/

      sqrt3=sqrt(three)

      do 1001 j=1,36
      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
      sss(i)=zero
 1005 if(efac(i).lt.expcut) sss(i)=zetafac(i)*base(i)*exp(-efac(i))

      nn=0

      do 1006 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
 1006 continue

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

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

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

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

 1010 continue

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

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

      do 1013 j=1,6
      do 1013 i=1,nn
 1013 c2(i,j)=zetarec(i)*dss(i,j)

      do 1014 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
 1014 gss(i,15)=gzaz(i)*fss(i,10)+c2(i,6)*three

      do 1015 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)
 1015 dps(i,18)=azbz(iicontr(i))*dss(i,6)+fss(i,10)

      do 1016 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)
 1016 fps(i,30)=azbz(iicontr(i))*fss(i,10)+gss(i,15)

      do 1017 i=1,nn
      dds(i,01)=axbx(iicontr(i))*dps(i,01)+fps(i,01)
      dds(i,02)=ayby(iicontr(i))*dps(i,01)+fps(i,04)
      dds(i,03)=azbz(iicontr(i))*dps(i,01)+fps(i,07)
      dds(i,04)=ayby(iicontr(i))*dps(i,02)+fps(i,05)
      dds(i,05)=azbz(iicontr(i))*dps(i,02)+fps(i,08)
      dds(i,06)=azbz(iicontr(i))*dps(i,03)+fps(i,09)
      dds(i,07)=axbx(iicontr(i))*dps(i,04)+fps(i,04)
      dds(i,08)=ayby(iicontr(i))*dps(i,04)+fps(i,10)
      dds(i,09)=azbz(iicontr(i))*dps(i,04)+fps(i,13)
      dds(i,10)=ayby(iicontr(i))*dps(i,05)+fps(i,11)
      dds(i,11)=azbz(iicontr(i))*dps(i,05)+fps(i,14)
      dds(i,12)=azbz(iicontr(i))*dps(i,06)+fps(i,15)
      dds(i,13)=axbx(iicontr(i))*dps(i,07)+fps(i,07)
      dds(i,14)=ayby(iicontr(i))*dps(i,07)+fps(i,13)
      dds(i,15)=azbz(iicontr(i))*dps(i,07)+fps(i,16)
      dds(i,16)=ayby(iicontr(i))*dps(i,08)+fps(i,14)
      dds(i,17)=azbz(iicontr(i))*dps(i,08)+fps(i,17)
      dds(i,18)=azbz(iicontr(i))*dps(i,09)+fps(i,18)
      dds(i,19)=axbx(iicontr(i))*dps(i,10)+fps(i,10)
      dds(i,20)=ayby(iicontr(i))*dps(i,10)+fps(i,19)
      dds(i,21)=azbz(iicontr(i))*dps(i,10)+fps(i,22)
      dds(i,22)=ayby(iicontr(i))*dps(i,11)+fps(i,20)
      dds(i,23)=azbz(iicontr(i))*dps(i,11)+fps(i,23)
      dds(i,24)=azbz(iicontr(i))*dps(i,12)+fps(i,24)
      dds(i,25)=axbx(iicontr(i))*dps(i,13)+fps(i,13)
      dds(i,26)=ayby(iicontr(i))*dps(i,13)+fps(i,22)
      dds(i,27)=azbz(iicontr(i))*dps(i,13)+fps(i,25)
      dds(i,28)=ayby(iicontr(i))*dps(i,14)+fps(i,23)
      dds(i,29)=azbz(iicontr(i))*dps(i,14)+fps(i,26)
      dds(i,30)=azbz(iicontr(i))*dps(i,15)+fps(i,27)
      dds(i,31)=axbx(iicontr(i))*dps(i,16)+fps(i,16)
      dds(i,32)=ayby(iicontr(i))*dps(i,16)+fps(i,25)
      dds(i,33)=azbz(iicontr(i))*dps(i,16)+fps(i,28)
      dds(i,34)=ayby(iicontr(i))*dps(i,17)+fps(i,26)
      dds(i,35)=azbz(iicontr(i))*dps(i,17)+fps(i,29)
 1017 dds(i,36)=azbz(iicontr(i))*dps(i,18)+fps(i,30)

      do 1018 j=1,36
cfpp$ nodepchk
      do 1018 i=1,nn
 1018 gather(iloc(i),j)=gather(iloc(i),j)+dds(i,j)*xcfits

 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 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 j=1,15
      do 2016 i=1,nn
 2016 c4(i,j)=zetarec(i)*gss(i,j)

      do 2017 i=1,nn
      qss(i,01)=gxax(i)*hss(i,01)+c4(i,01)*five
      qss(i,02)=gyay(i)*hss(i,01)
      qss(i,03)=gzaz(i)*hss(i,01)
      qss(i,04)=gyay(i)*hss(i,02)+c4(i,01)
      qss(i,05)=gzaz(i)*hss(i,02)
      qss(i,06)=gzaz(i)*hss(i,03)+c4(i,01)
      qss(i,07)=gyay(i)*hss(i,04)+c4(i,02)*two
      qss(i,08)=gzaz(i)*hss(i,04)
      qss(i,09)=gzaz(i)*hss(i,05)+c4(i,02)
      qss(i,10)=gzaz(i)*hss(i,06)+c4(i,03)*two
      qss(i,11)=gyay(i)*hss(i,07)+c4(i,04)*three
      qss(i,12)=gzaz(i)*hss(i,07)
      qss(i,13)=gzaz(i)*hss(i,08)+c4(i,04)
      qss(i,14)=gzaz(i)*hss(i,09)+c4(i,05)*two
      qss(i,15)=gzaz(i)*hss(i,10)+c4(i,06)*three
      qss(i,16)=gyay(i)*hss(i,11)+c4(i,07)*four
      qss(i,17)=gzaz(i)*hss(i,11)
      qss(i,18)=gzaz(i)*hss(i,12)+c4(i,07)
      qss(i,19)=gzaz(i)*hss(i,13)+c4(i,08)*two
      qss(i,20)=gzaz(i)*hss(i,14)+c4(i,09)*three
      qss(i,21)=gzaz(i)*hss(i,15)+c4(i,10)*four
      qss(i,22)=gyay(i)*hss(i,16)+c4(i,11)*five
      qss(i,23)=gzaz(i)*hss(i,16)
      qss(i,24)=gzaz(i)*hss(i,17)+c4(i,11)
      qss(i,25)=gzaz(i)*hss(i,18)+c4(i,12)*two
      qss(i,26)=gzaz(i)*hss(i,19)+c4(i,13)*three
      qss(i,27)=gzaz(i)*hss(i,20)+c4(i,14)*four
 2017 qss(i,28)=gzaz(i)*hss(i,21)+c4(i,15)*five

      call xcordd1(m,n,nn,iicontr,axbx,ayby,azbz,axcx,aycy,azcz,dss,
     &             fss,gss,hss,qss,dps,dds,ddp,ddd,fps,fds,fdp,gps,gds,
     &             hps)

      do 2018 j=1,36
cfpp$ nodepchk
      do 2018 i=1,nn
 2018 gather(iloc(i),j)=gather(iloc(i),j)
     &                 +dds(i,j)*xcfits
     &                 +ddp(i,3*(j-1)+1)*xcfitp(1)
     &                 +ddp(i,3*(j-1)+2)*xcfitp(2)
     &                 +ddp(i,3*(j-1)+3)*xcfitp(3)
     &                 +ddd(i,6*(j-1)+1)*xcfitd(1)
     &                 +ddd(i,6*(j-1)+2)*xcfitd(2)
     &                 +ddd(i,6*(j-1)+3)*xcfitd(3)
     &                 +ddd(i,6*(j-1)+4)*xcfitd(4)
     &                 +ddd(i,6*(j-1)+5)*xcfitd(5)
     &                 +ddd(i,6*(j-1)+6)*xcfitd(6)

 2001 continue

      index(01)=01
      index(02)=07
      index(03)=08
      index(04)=13
      index(05)=14
      index(06)=15
      index(07)=19
      index(08)=20
      index(09)=21
      index(10)=22
      index(11)=25
      index(12)=26
      index(13)=27
      index(14)=28
      index(15)=29
      index(16)=31
      index(17)=32
      index(18)=33
      index(19)=34
      index(20)=35
      index(21)=36

      do 1019 i=1,n
cfpp$ nodepchk
      do 1019 j=1,21
      jj=index(j)
 1019 fock3(inddyy(icontr(i),jj))=fock3(inddyy(icontr(i),jj))
     &                           +gather(i,jj)

      index(01)=02
      index(02)=03
      index(03)=04
      index(04)=05
      index(05)=06
      index(06)=09
      index(07)=10
      index(08)=11
      index(09)=12
      index(10)=16
      index(11)=17
      index(12)=18
      index(13)=23
      index(14)=24
      index(15)=30

      do 1020 i=1,n
      if(iminj(icontr(i)).ne.0) then

cfpp$ nodepchk
      do 1021 j=1,15
      jj=index(j)
 1021 fock3(inddyy(icontr(i),jj))=fock3(inddyy(icontr(i),jj))
     &                           +gather(i,jj)

                                endif
 1020 continue

      return
      end
