      subroutine boys_mm(nreal,ninteger,ncontractions,nconts,ncontp,
     &                   ncontd,nshels,nshelp,nsheld,mtloca,iprimvc,
     &                   nprimvc,icentvc,iwkvec,coord,alpha,coeff,
     &                   x_moment,y_moment,z_moment,wkvec)

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 nshels(*),nshelp(*),nsheld(*),mtloca(*),iprimvc(*),
     &          nprimvc(*),icentvc(*),iwkvec(*)

      dimension coord(3,*),alpha(*),coeff(*),x_moment(ncontractions,*),
     &          y_moment(ncontractions,*),z_moment(ncontractions,*),
     &          wkvec(*)

      data zero,cutoff/0.0d0,1.0d+35/

      do 1001 j=1,ncontractions
      do 1001 i=1,ncontractions
      x_moment(i,j)=zero
      y_moment(i,j)=zero
 1001 z_moment(i,j)=zero

      ncposs=(nconts*(nconts+1))/2
      ncpopp=(ncontp*(ncontp+1))/2
      ncpodd=(ncontd*(ncontd+1))/2

      ncpops=ncontp*nconts
      ncpods=ncontd*nconts
      ncpodp=ncontd*ncontp

      npposs=0
      nppops=0
      nppopp=0
      nppods=0
      nppodp=0
      nppodd=0

      do 1002 i=1,nconts
      do 1002 j=1,i
 1002 npposs=npposs+nprimvc(nshels(i))*nprimvc(nshels(j))

      do 1003 i=1,ncontp
      do 1003 j=1,nconts
 1003 nppops=nppops+nprimvc(nshelp(i))*nprimvc(nshels(j))

      do 1004 i=1,ncontp
      do 1004 j=1,i
 1004 nppopp=nppopp+nprimvc(nshelp(i))*nprimvc(nshelp(j))

      do 1005 i=1,ncontd
      do 1005 j=1,nconts
 1005 nppods=nppods+nprimvc(nsheld(i))*nprimvc(nshels(j))

      do 1006 i=1,ncontd
      do 1006 j=1,ncontp
 1006 nppodp=nppodp+nprimvc(nsheld(i))*nprimvc(nshelp(j))

      do 1007 i=1,ncontd
      do 1007 j=1,i
 1007 nppodd=nppodd+nprimvc(nsheld(i))*nprimvc(nsheld(j))

      ja=1
      jb=ja+npposs
      jc=jb+ncposs

      if(jc.gt.ninteger) stop 'boys_mm ss integer overflow'

      je=1
      jf=je+npposs
      jg=jf+npposs
      jh=jg+npposs
      ji=jh+npposs
      jj=ji+npposs
      jk=jj+npposs
      jl=jk+npposs
      jm=jl+npposs
      jn=jm+npposs
      jo=jn+npposs
      jp=jo+npposs
      jq=jp+npposs

      if(jq.gt.nreal) stop 'boys_mm ss real overflow'

      call setupss(ncposs,n,nconts,nshels,mtloca,iprimvc,nprimvc,
     &             icentvc,iwkvec(ja),iwkvec(jb),cutoff,coord,alpha,
     &             coeff,wkvec(je),wkvec(jf),wkvec(jg),wkvec(jh),
     &             wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),wkvec(jm),
     &             wkvec(jn),wkvec(jo),wkvec(jp))

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n

      ka=jq+1
      kb=ka+n*3

      if(kb.gt.nreal) stop 'boys_mm ss real overflow'

      call boys_ss(n,iwkvec(ja),iwkvec(jb),wkvec(je),wkvec(jn),
     &             wkvec(jo),wkvec(jp),wkvec(jm),wkvec(ka),x_moment,
     &             y_moment,z_moment)

      if(ncontp.eq.0) goto 1008

      ja=1
      jb=ja+nppops
      jc=jb+ncpops*3

      if(jc.gt.ninteger) stop 'boys_mm ps integer overflow'

      je=1
      jf=je+nppops
      jg=jf+nppops
      jh=jg+nppops
      ji=jh+nppops
      jj=ji+nppops
      jk=jj+nppops
      jl=jk+nppops
      jm=jl+nppops
      jn=jm+nppops
      jo=jn+nppops
      jp=jo+nppops
      jq=jp+nppops
      jr=jq+nppops
      js=jr+nppops
      jt=js+nppops

      if(jt.gt.nreal) stop 'boys_mm ps real overflow'

      call setupps(ncpops,n,ncontp,nconts,nshelp,nshels,mtloca,iprimvc,
     &             nprimvc,icentvc,iwkvec(ja),iwkvec(jb),cutoff,coord,
     &             alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),wkvec(jh),
     &             wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),wkvec(jm),
     &             wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),wkvec(jr),
     &             wkvec(js))

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))
      call stringcopy(n,wkvec(12*n+1),wkvec(jq))
      call stringcopy(n,wkvec(13*n+1),wkvec(jr))
      call stringcopy(n,wkvec(14*n+1),wkvec(js))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n

      ka=jt+1
      kb=ka+n*3
      kc=kb+n*9
      kd=kc+n

      if(kd.gt.nreal) stop 'boys_mm ps real overflow'

      call boys_ps(ncpops,n,iwkvec(ja),iwkvec(jb),wkvec(je),wkvec(jn),
     &             wkvec(jo),wkvec(jp),wkvec(jq),wkvec(jr),wkvec(js),
     &             wkvec(jm),wkvec(ka),wkvec(kb),wkvec(kc),x_moment,
     &             y_moment,z_moment)

      ja=1
      jb=ja+ncpopp
      jc=jb+nppopp
      jd=jc+ncpopp*9

      if(jd.gt.ninteger) stop 'boys_mm pp integer overflow'

      je=1
      jf=je+nppopp
      jg=jf+nppopp
      jh=jg+nppopp
      ji=jh+nppopp
      jj=ji+nppopp
      jk=jj+nppopp
      jl=jk+nppopp
      jm=jl+nppopp
      jn=jm+nppopp
      jo=jn+nppopp
      jp=jo+nppopp
      jq=jp+nppopp
      jr=jq+nppopp
      js=jr+nppopp
      jt=js+nppopp
      ju=jt+nppopp
      jv=ju+nppopp
      jw=jv+nppopp
      jx=jw+ncpopp
      jy=jx+ncpopp
      jz=jy+ncpopp

      if(jz.gt.nreal) stop 'boys_mm pp real overflow'

      call setuppp(ncpopp,n,ncontp,nshelp,mtloca,iprimvc,nprimvc,
     &             icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),cutoff,
     &             coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &             wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &             wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &             wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &             wkvec(jw),wkvec(jx),wkvec(jy))

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))
      call stringcopy(n,wkvec(12*n+1),wkvec(jq))
      call stringcopy(n,wkvec(13*n+1),wkvec(jr))
      call stringcopy(n,wkvec(14*n+1),wkvec(js))
      call stringcopy(n,wkvec(15*n+1),wkvec(jt))
      call stringcopy(n,wkvec(16*n+1),wkvec(ju))
      call stringcopy(n,wkvec(17*n+1),wkvec(jv))
      call stringcopy(ncpopp,wkvec(18*n+0*ncpopp+1),wkvec(jw))
      call stringcopy(ncpopp,wkvec(18*n+1*ncpopp+1),wkvec(jx))
      call stringcopy(ncpopp,wkvec(18*n+2*ncpopp+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpopp
      jy=jx+ncpopp
      jz=jy+ncpopp

      ka=jz+1
      kb=ka+n*3
      kc=kb+n*9
      kd=kc+n*27
      ke=kd+n*3
      kf=ke+n
      kg=kf+n*3
      kh=kg+n*3

      if(kh.gt.nreal) stop 'boys_mm pp real overflow'

      call boys_pp(ncpopp,n,iwkvec(ja),iwkvec(jb),iwkvec(jc),wkvec(je),
     &             wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),wkvec(jr),
     &             wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jm),
     &             wkvec(ka),wkvec(kb),wkvec(kc),wkvec(kd),wkvec(ke),
     &             wkvec(kf),wkvec(kg),x_moment,y_moment,z_moment)

      if(ncontd.eq.0) goto 1008

      ja=1
      jb=ja+nppods
      jc=jb+ncpods*6

      if(jc.gt.ninteger) stop 'boys_mm ds integer overflow'

      je=1
      jf=je+nppods
      jg=jf+nppods
      jh=jg+nppods
      ji=jh+nppods
      jj=ji+nppods
      jk=jj+nppods
      jl=jk+nppods
      jm=jl+nppods
      jn=jm+nppods
      jo=jn+nppods
      jp=jo+nppods
      jq=jp+nppods
      jr=jq+nppods
      js=jr+nppods
      jt=js+nppods

      if(jt.gt.nreal) stop 'boys_mm ds real overflow'

      call setupds(ncpods,n,ncontd,nconts,nsheld,nshels,mtloca,iprimvc,
     &             nprimvc,icentvc,iwkvec(ja),iwkvec(jb),cutoff,coord,
     &             alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),wkvec(jh),
     &             wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),wkvec(jm),
     &             wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),wkvec(jr),
     &             wkvec(js))

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))
      call stringcopy(n,wkvec(12*n+1),wkvec(jq))
      call stringcopy(n,wkvec(13*n+1),wkvec(jr))
      call stringcopy(n,wkvec(14*n+1),wkvec(js))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n

      ka=jt+1
      kb=ka+n*3
      kc=kb+n*9
      kd=kc+n*18
      ke=kd+n*3
      kf=ke+n
      kg=kf+n*3
      kh=kg+n*3

      if(kh.gt.nreal) stop 'boys_mm ds real overflow'

      call boys_ds(ncpods,n,iwkvec(ja),iwkvec(jb),wkvec(je),wkvec(jn),
     &             wkvec(jo),wkvec(jp),wkvec(jq),wkvec(jr),wkvec(js),
     &             wkvec(jm),wkvec(ka),wkvec(kb),wkvec(kc),wkvec(kd),
     &             wkvec(ke),wkvec(kf),wkvec(kg),x_moment,y_moment,
     &             z_moment)

      ja=1
      jb=ja+nppodp
      jc=jb+ncpodp*18

      if(jc.gt.ninteger) stop 'boys_mm dp integer overflow'

      je=1
      jf=je+nppodp
      jg=jf+nppodp
      jh=jg+nppodp
      ji=jh+nppodp
      jj=ji+nppodp
      jk=jj+nppodp
      jl=jk+nppodp
      jm=jl+nppodp
      jn=jm+nppodp
      jo=jn+nppodp
      jp=jo+nppodp
      jq=jp+nppodp
      jr=jq+nppodp
      js=jr+nppodp
      jt=js+nppodp
      ju=jt+nppodp
      jv=ju+nppodp
      jw=jv+nppodp
      jx=jw+ncpodp
      jy=jx+ncpodp
      jz=jy+ncpodp

      if(jz.gt.nreal) stop 'boys_mm dp real overflow'

      call setupdp(ncpodp,n,ncontd,ncontp,nsheld,nshelp,mtloca,iprimvc,
     &             nprimvc,icentvc,iwkvec(ja),iwkvec(jb),cutoff,coord,
     &             alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),wkvec(jh),
     &             wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),wkvec(jm),
     &             wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),wkvec(jr),
     &             wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jw),
     &             wkvec(jx),wkvec(jy))

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))
      call stringcopy(n,wkvec(12*n+1),wkvec(jq))
      call stringcopy(n,wkvec(13*n+1),wkvec(jr))
      call stringcopy(n,wkvec(14*n+1),wkvec(js))
      call stringcopy(n,wkvec(15*n+1),wkvec(jt))
      call stringcopy(n,wkvec(16*n+1),wkvec(ju))
      call stringcopy(n,wkvec(17*n+1),wkvec(jv))
      call stringcopy(ncpodp,wkvec(18*n+0*ncpodp+1),wkvec(jw))
      call stringcopy(ncpodp,wkvec(18*n+1*ncpodp+1),wkvec(jx))
      call stringcopy(ncpodp,wkvec(18*n+2*ncpodp+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpodp
      jy=jx+ncpodp
      jz=jy+ncpodp

      ka=jz+1
      kb=ka+n*3
      kc=kb+n*9
      kd=kc+n*18
      ke=kd+n*54
      kf=ke+n*3
      kg=kf+n*6
      kh=kg+n
      ki=kh+n*6
      kj=ki+n*9

      if(kj.gt.nreal) stop 'boys_mm dp real overflow'

      call boys_dp(ncpodp,n,iwkvec(ja),iwkvec(jb),wkvec(je),wkvec(jn),
     &             wkvec(jo),wkvec(jp),wkvec(jq),wkvec(jr),wkvec(js),
     &             wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jm),wkvec(ka),
     &             wkvec(kb),wkvec(kc),wkvec(kd),wkvec(ke),wkvec(kf),
     &             wkvec(kg),wkvec(kh),wkvec(ki),x_moment,y_moment,
     &             z_moment)

      ja=1
      jb=ja+ncpodd
      jc=jb+nppodd
      jd=jc+ncpodd*36

      if(jd.gt.ninteger) stop 'boys_mm dd integer overflow'

      je=1
      jf=je+nppodd
      jg=jf+nppodd
      jh=jg+nppodd
      ji=jh+nppodd
      jj=ji+nppodd
      jk=jj+nppodd
      jl=jk+nppodd
      jm=jl+nppodd
      jn=jm+nppodd
      jo=jn+nppodd
      jp=jo+nppodd
      jq=jp+nppodd
      jr=jq+nppodd
      js=jr+nppodd
      jt=js+nppodd
      ju=jt+nppodd
      jv=ju+nppodd
      jw=jv+nppodd
      jx=jw+ncpodd
      jy=jx+ncpodd
      jz=jy+ncpodd

      if(jz.gt.nreal) stop 'boys_mm dd real overflow'

      call setupdd(ncpodd,n,ncontd,nsheld,mtloca,iprimvc,nprimvc,
     &             icentvc,iwkvec(ja),iwkvec(jb),iwkvec(jc),cutoff,
     &             coord,alpha,coeff,wkvec(je),wkvec(jf),wkvec(jg),
     &             wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),wkvec(jl),
     &             wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),
     &             wkvec(jr),wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),
     &             wkvec(jw),wkvec(jx),wkvec(jy))

      call stringcopy(n,wkvec(00*n+1),wkvec(je))
      call stringcopy(n,wkvec(01*n+1),wkvec(jf))
      call stringcopy(n,wkvec(02*n+1),wkvec(jg))
      call stringcopy(n,wkvec(03*n+1),wkvec(jh))
      call stringcopy(n,wkvec(04*n+1),wkvec(ji))
      call stringcopy(n,wkvec(05*n+1),wkvec(jj))
      call stringcopy(n,wkvec(06*n+1),wkvec(jk))
      call stringcopy(n,wkvec(07*n+1),wkvec(jl))
      call stringcopy(n,wkvec(08*n+1),wkvec(jm))
      call stringcopy(n,wkvec(09*n+1),wkvec(jn))
      call stringcopy(n,wkvec(10*n+1),wkvec(jo))
      call stringcopy(n,wkvec(11*n+1),wkvec(jp))
      call stringcopy(n,wkvec(12*n+1),wkvec(jq))
      call stringcopy(n,wkvec(13*n+1),wkvec(jr))
      call stringcopy(n,wkvec(14*n+1),wkvec(js))
      call stringcopy(n,wkvec(15*n+1),wkvec(jt))
      call stringcopy(n,wkvec(16*n+1),wkvec(ju))
      call stringcopy(n,wkvec(17*n+1),wkvec(jv))
      call stringcopy(ncpodd,wkvec(18*n+0*ncpodd+1),wkvec(jw))
      call stringcopy(ncpodd,wkvec(18*n+1*ncpodd+1),wkvec(jx))
      call stringcopy(ncpodd,wkvec(18*n+2*ncpodd+1),wkvec(jy))

      je=1
      jf=je+n
      jg=jf+n
      jh=jg+n
      ji=jh+n
      jj=ji+n
      jk=jj+n
      jl=jk+n
      jm=jl+n
      jn=jm+n
      jo=jn+n
      jp=jo+n
      jq=jp+n
      jr=jq+n
      js=jr+n
      jt=js+n
      ju=jt+n
      jv=ju+n
      jw=jv+n
      jx=jw+ncpodd
      jy=jx+ncpodd
      jz=jy+ncpodd

      ka=jz+1
      kb=ka+n*3
      kc=kb+n*9
      kd=kc+n*18
      ke=kd+n*27
      kf=ke+n*54
      kg=kf+n*108
      kh=kg+n*3
      ki=kh+n*6
      kj=ki+n*18
      kk=kj+n
      kl=kk+n*18
      km=kl+n*18
      kn=km+n*27

      if(kn.gt.nreal) stop 'boys_mm dd real overflow'

      call boys_dd(ncpodd,n,iwkvec(ja),iwkvec(jb),iwkvec(jc),wkvec(je),
     &             wkvec(jn),wkvec(jo),wkvec(jp),wkvec(jq),wkvec(jr),
     &             wkvec(js),wkvec(jt),wkvec(ju),wkvec(jv),wkvec(jm),
     &             wkvec(ka),wkvec(kb),wkvec(kc),wkvec(kd),wkvec(ke),
     &             wkvec(kf),wkvec(kg),wkvec(kh),wkvec(ki),wkvec(kj),
     &             wkvec(kk),wkvec(kl),wkvec(km),x_moment,y_moment,
     &             z_moment)

 1008 continue

      ndim=(ncontractions*(ncontractions+1))/2

      call stringcopy(ndim,wkvec,x_moment)

      ij=0
      do 1009 i=1,ncontractions
      do 1009 j=1,i
      ij=ij+1
      x_moment(i,j)=wkvec(ij)
 1009 x_moment(j,i)=wkvec(ij)

      call stringcopy(ndim,wkvec,y_moment)

      ij=0
      do 1010 i=1,ncontractions
      do 1010 j=1,i
      ij=ij+1
      y_moment(i,j)=wkvec(ij)
 1010 y_moment(j,i)=wkvec(ij)

      call stringcopy(ndim,wkvec,z_moment)

      ij=0
      do 1011 i=1,ncontractions
      do 1011 j=1,i
      ij=ij+1
      z_moment(i,j)=wkvec(ij)
 1011 z_moment(j,i)=wkvec(ij)

      return
      end
