      subroutine compg2(iatom,ncntrt,npts,nd,natoms,nconts,ncontp,
     &                  ncontd,ilfunc,icfunc,ngauss,nptrs,nptrp,nptrd,
     &                  radius,cent,alpha,coeff,x,y,z,qf,qx,qy,qz,q1,q2,
     &                  q3,q4,q5,q6,extent1,extent3,dcc,xx,yy,zz,rr,
     &                  rralpha,s,px,py,pz,dxx,dxy,dxz,dyy,dyz,dzz,fxxx,
     &                  fxxy,fxxz,fxyy,fxyz,fxzz,fyyy,fyyz,fyzz,fzzz,
     &                  gxxxx,gxxxy,gxxxz,gxxyy,gxxyz,gxxzz,gxyyy,gxyyz,
     &                  gxyzz,gxzzz,gyyyy,gyyyz,gyyzz,gyzzz,gzzzz)

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 ilfunc(*),icfunc(*),ngauss(*),nptrs(*),nptrp(*),nptrd(*)

      dimension cent(3,*),alpha(*),coeff(*)

      dimension x(*),y(*),z(*),qf(nd,*),qx(nd,*),qy(nd,*),qz(nd,*)
      dimension q1(nd,*),q2(nd,*),q3(nd,*),q4(nd,*),q5(nd,*),q6(nd,*)

      dimension extent1(*),extent3(*),dcc(natoms,*)
      dimension xx(*),yy(*),zz(*),rr(*),rralpha(*)
      dimension s(*),px(*),py(*),pz(*)
      dimension dxx(*),dxy(*),dxz(*),dyy(*),dyz(*),dzz(*)
      dimension fxxx(*),fxxy(*),fxxz(*),fxyy(*),fxyz(*),
     &          fxzz(*),fyyy(*),fyyz(*),fyzz(*),fzzz(*)
      dimension gxxxx(*),gxxxy(*),gxxxz(*),gxxyy(*),gxxyz(*),
     &          gxxzz(*),gxyyy(*),gxyyz(*),gxyzz(*),gxzzz(*),
     &          gyyyy(*),gyyyz(*),gyyzz(*),gyzzz(*),gzzzz(*)

      data zero,one,two,three,five,expcut/0.0d0,1.0d0,2.0d0,
     &                                    3.0d0,5.0d0,50.0d0/

      rsqrt3=one/sqrt(three)

      do 1001 i=1,ncntrt
 1001 extent3(i)=zero

      k=0

      if(nconts.eq.0) return

      do 1002 i=1,nconts

      ilabel=nptrs(i)

      nc=icfunc(ilabel)

      distmin=dcc(nc,iatom)-radius
      if(distmin.gt.extent1(i)) goto 1002

      iptr1=ilfunc(ilabel)
      iptr2=iptr1+ngauss(ilabel)-1

      do 1003 j=1,npts
      qf(j,k+1)=zero
      qx(j,k+1)=zero
      qy(j,k+1)=zero
      qz(j,k+1)=zero
      q1(j,k+1)=zero
      q2(j,k+1)=zero
      q3(j,k+1)=zero
      q4(j,k+1)=zero
      q5(j,k+1)=zero
 1003 q6(j,k+1)=zero

      do 1004 j=1,npts
      xx(j)=x(j)-cent(1,nc)
      yy(j)=y(j)-cent(2,nc)
 1004 zz(j)=z(j)-cent(3,nc)

      do 1005 j=1,npts
 1005 rr(j)=xx(j)*xx(j)+yy(j)*yy(j)+zz(j)*zz(j)

      do 1006 nprimitive=iptr1,iptr2

      twozeta=two*alpha(nprimitive)

      do 1007 j=1,npts
 1007 rralpha(j)=rr(j)*alpha(nprimitive)

      do 1008 j=1,npts
      if(rralpha(j).lt.expcut) then

                       s(j)=coeff(nprimitive)*exp(-rralpha(j))

                       qf(j,k+1)=qf(j,k+1)+s(j)

                       s(j)=s(j)*twozeta

                       px(j)=xx(j)*s(j)
                       py(j)=yy(j)*s(j)
                       pz(j)=zz(j)*s(j)

                       qx(j,k+1)=qx(j,k+1)-px(j)
                       qy(j,k+1)=qy(j,k+1)-py(j)
                       qz(j,k+1)=qz(j,k+1)-pz(j)

                       px(j)=px(j)*twozeta
                       py(j)=py(j)*twozeta
                       pz(j)=pz(j)*twozeta

                       q1(j,k+1)=q1(j,k+1)+xx(j)*px(j)-s(j)
                       q2(j,k+1)=q2(j,k+1)+xx(j)*py(j)
                       q3(j,k+1)=q3(j,k+1)+xx(j)*pz(j)
                       q4(j,k+1)=q4(j,k+1)+yy(j)*py(j)-s(j)
                       q5(j,k+1)=q5(j,k+1)+yy(j)*pz(j)
                       q6(j,k+1)=q6(j,k+1)+zz(j)*pz(j)-s(j)

                               endif
 1008 continue

 1006 continue

      extent3(k+1)=sumsum(npts,qf(1,k+1))/real64(npts)

 1002 k=k+1

      if(ncontp.eq.0) return

      do 1009 i=1,ncontp

      ilabel=nptrp(i)

      nc=icfunc(ilabel)

      distmin=dcc(nc,iatom)-radius
      if(distmin.gt.extent1(nconts+i)) goto 1009

      iptr1=ilfunc(ilabel)
      iptr2=iptr1+ngauss(ilabel)-1

      do 1010 j=1,npts
      qf(j,k+1)=zero
      qx(j,k+1)=zero
      qy(j,k+1)=zero
      qz(j,k+1)=zero
      q1(j,k+1)=zero
      q2(j,k+1)=zero
      q3(j,k+1)=zero
      q4(j,k+1)=zero
      q5(j,k+1)=zero
      q6(j,k+1)=zero
      qf(j,k+2)=zero
      qx(j,k+2)=zero
      qy(j,k+2)=zero
      qz(j,k+2)=zero
      q1(j,k+2)=zero
      q2(j,k+2)=zero
      q3(j,k+2)=zero
      q4(j,k+2)=zero
      q5(j,k+2)=zero
      q6(j,k+2)=zero
      qf(j,k+3)=zero
      qx(j,k+3)=zero
      qy(j,k+3)=zero
      qz(j,k+3)=zero
      q1(j,k+3)=zero
      q2(j,k+3)=zero
      q3(j,k+3)=zero
      q4(j,k+3)=zero
      q5(j,k+3)=zero
 1010 q6(j,k+3)=zero

      do 1011 j=1,npts
      xx(j)=x(j)-cent(1,nc)
      yy(j)=y(j)-cent(2,nc)
 1011 zz(j)=z(j)-cent(3,nc)

      do 1012 j=1,npts
 1012 rr(j)=xx(j)*xx(j)+yy(j)*yy(j)+zz(j)*zz(j)

      do 1013 nprimitive=iptr1,iptr2

      twozeta=two*alpha(nprimitive)

      do 1014 j=1,npts
 1014 rralpha(j)=rr(j)*alpha(nprimitive)

      do 1015 j=1,npts
      if(rralpha(j).lt.expcut) then

                       s(j)=coeff(nprimitive)*exp(-rralpha(j))

                       px(j)=xx(j)*s(j)
                       py(j)=yy(j)*s(j)
                       pz(j)=zz(j)*s(j)

                       qf(j,k+1)=qf(j,k+1)+px(j)
                       qf(j,k+2)=qf(j,k+2)+py(j)
                       qf(j,k+3)=qf(j,k+3)+pz(j)

                       px(j)=px(j)*twozeta
                       py(j)=py(j)*twozeta
                       pz(j)=pz(j)*twozeta

                       dxx(j)=xx(j)*px(j)
                       dxy(j)=xx(j)*py(j)
                       dxz(j)=xx(j)*pz(j)
                       dyy(j)=yy(j)*py(j)
                       dyz(j)=yy(j)*pz(j)
                       dzz(j)=zz(j)*pz(j)

                       qx(j,k+1)=qx(j,k+1)-dxx(j)+s(j)
                       qy(j,k+1)=qy(j,k+1)-dxy(j)
                       qz(j,k+1)=qz(j,k+1)-dxz(j)

                       qx(j,k+2)=qx(j,k+2)-dxy(j)
                       qy(j,k+2)=qy(j,k+2)-dyy(j)+s(j)
                       qz(j,k+2)=qz(j,k+2)-dyz(j)

                       qx(j,k+3)=qx(j,k+3)-dxz(j)
                       qy(j,k+3)=qy(j,k+3)-dyz(j)
                       qz(j,k+3)=qz(j,k+3)-dzz(j)+s(j)

                       dxx(j)=dxx(j)*twozeta
                       dxy(j)=dxy(j)*twozeta
                       dxz(j)=dxz(j)*twozeta
                       dyy(j)=dyy(j)*twozeta
                       dyz(j)=dyz(j)*twozeta
                       dzz(j)=dzz(j)*twozeta

                       fxxx(j)=xx(j)*dxx(j)
                       fxxy(j)=xx(j)*dxy(j)
                       fxxz(j)=xx(j)*dxz(j)
                       fxyy(j)=xx(j)*dyy(j)
                       fxyz(j)=xx(j)*dyz(j)
                       fxzz(j)=xx(j)*dzz(j)
                       fyyy(j)=yy(j)*dyy(j)
                       fyyz(j)=yy(j)*dyz(j)
                       fyzz(j)=yy(j)*dzz(j)
                       fzzz(j)=zz(j)*dzz(j)

                       q1(j,k+1)=q1(j,k+1)+fxxx(j)-px(j)-px(j)-px(j)
                       q2(j,k+1)=q2(j,k+1)+fxxy(j)-py(j)
                       q3(j,k+1)=q3(j,k+1)+fxxz(j)-pz(j)
                       q4(j,k+1)=q4(j,k+1)+fxyy(j)-px(j)
                       q5(j,k+1)=q5(j,k+1)+fxyz(j)
                       q6(j,k+1)=q6(j,k+1)+fxzz(j)-px(j)

                       q1(j,k+2)=q1(j,k+2)+fxxy(j)-py(j)
                       q2(j,k+2)=q2(j,k+2)+fxyy(j)-px(j)
                       q3(j,k+2)=q3(j,k+2)+fxyz(j)
                       q4(j,k+2)=q4(j,k+2)+fyyy(j)-py(j)-py(j)-py(j)
                       q5(j,k+2)=q5(j,k+2)+fyyz(j)-pz(j)
                       q6(j,k+2)=q6(j,k+2)+fyzz(j)-py(j)

                       q1(j,k+3)=q1(j,k+3)+fxxz(j)-pz(j)
                       q2(j,k+3)=q2(j,k+3)+fxyz(j)
                       q3(j,k+3)=q3(j,k+3)+fxzz(j)-px(j)
                       q4(j,k+3)=q4(j,k+3)+fyyz(j)-pz(j)
                       q5(j,k+3)=q5(j,k+3)+fyzz(j)-py(j)
                       q6(j,k+3)=q6(j,k+3)+fzzz(j)-pz(j)-pz(j)-pz(j)

                               endif
 1015 continue

 1013 continue

      extent3(k+1)=sumsum(npts,qf(1,k+1))/real64(npts)
      extent3(k+2)=sumsum(npts,qf(1,k+2))/real64(npts)
      extent3(k+3)=sumsum(npts,qf(1,k+3))/real64(npts)

 1009 k=k+3

      if(ncontd.eq.0) return

      do 1016 i=1,ncontd

      ilabel=nptrd(i)

      nc=icfunc(ilabel)

      distmin=dcc(nc,iatom)-radius
      if(distmin.gt.extent1(nconts+ncontp+i)) goto 1016

      iptr1=ilfunc(ilabel)
      iptr2=iptr1+ngauss(ilabel)-1

      do 1017 j=1,npts
      qf(j,k+1)=zero
      qx(j,k+1)=zero
      qy(j,k+1)=zero
      qz(j,k+1)=zero
      q1(j,k+1)=zero
      q2(j,k+1)=zero
      q3(j,k+1)=zero
      q4(j,k+1)=zero
      q5(j,k+1)=zero
      q6(j,k+1)=zero
      qf(j,k+2)=zero
      qx(j,k+2)=zero
      qy(j,k+2)=zero
      qz(j,k+2)=zero
      q1(j,k+2)=zero
      q2(j,k+2)=zero
      q3(j,k+2)=zero
      q4(j,k+2)=zero
      q5(j,k+2)=zero
      q6(j,k+2)=zero
      qf(j,k+3)=zero
      qx(j,k+3)=zero
      qy(j,k+3)=zero
      qz(j,k+3)=zero
      q1(j,k+3)=zero
      q2(j,k+3)=zero
      q3(j,k+3)=zero
      q4(j,k+3)=zero
      q5(j,k+3)=zero
      q6(j,k+3)=zero
      qf(j,k+4)=zero
      qx(j,k+4)=zero
      qy(j,k+4)=zero
      qz(j,k+4)=zero
      q1(j,k+4)=zero
      q2(j,k+4)=zero
      q3(j,k+4)=zero
      q4(j,k+4)=zero
      q5(j,k+4)=zero
      q6(j,k+4)=zero
      qf(j,k+5)=zero
      qx(j,k+5)=zero
      qy(j,k+5)=zero
      qz(j,k+5)=zero
      q1(j,k+5)=zero
      q2(j,k+5)=zero
      q3(j,k+5)=zero
      q4(j,k+5)=zero
      q5(j,k+5)=zero
      q6(j,k+5)=zero
      qf(j,k+6)=zero
      qx(j,k+6)=zero
      qy(j,k+6)=zero
      qz(j,k+6)=zero
      q1(j,k+6)=zero
      q2(j,k+6)=zero
      q3(j,k+6)=zero
      q4(j,k+6)=zero
      q5(j,k+6)=zero
 1017 q6(j,k+6)=zero

      do 1018 j=1,npts
      xx(j)=x(j)-cent(1,nc)
      yy(j)=y(j)-cent(2,nc)
 1018 zz(j)=z(j)-cent(3,nc)

      do 1019 j=1,npts
 1019 rr(j)=xx(j)*xx(j)+yy(j)*yy(j)+zz(j)*zz(j)

      do 1020 nprimitive=iptr1,iptr2

      twozeta=two*alpha(nprimitive)

      do 1021 j=1,npts
 1021 rralpha(j)=rr(j)*alpha(nprimitive)

      do 1022 j=1,npts
      if(rralpha(j).lt.expcut) then

            s(j)=coeff(nprimitive)*exp(-rralpha(j))

            px(j)=xx(j)*s(j)
            py(j)=yy(j)*s(j)
            pz(j)=zz(j)*s(j)

            dxx(j)=xx(j)*px(j)
            dxy(j)=xx(j)*py(j)
            dxz(j)=xx(j)*pz(j)
            dyy(j)=yy(j)*py(j)
            dyz(j)=yy(j)*pz(j)
            dzz(j)=zz(j)*pz(j)

            qf(j,k+1)=qf(j,k+1)+dxx(j)
            qf(j,k+2)=qf(j,k+2)+dxy(j)
            qf(j,k+3)=qf(j,k+3)+dxz(j)
            qf(j,k+4)=qf(j,k+4)+dyy(j)
            qf(j,k+5)=qf(j,k+5)+dyz(j)
            qf(j,k+6)=qf(j,k+6)+dzz(j)

            dxx(j)=dxx(j)*twozeta
            dxy(j)=dxy(j)*twozeta
            dxz(j)=dxz(j)*twozeta
            dyy(j)=dyy(j)*twozeta
            dyz(j)=dyz(j)*twozeta
            dzz(j)=dzz(j)*twozeta

            fxxx(j)=xx(j)*dxx(j)
            fxxy(j)=xx(j)*dxy(j)
            fxxz(j)=xx(j)*dxz(j)
            fxyy(j)=xx(j)*dyy(j)
            fxyz(j)=xx(j)*dyz(j)
            fxzz(j)=xx(j)*dzz(j)
            fyyy(j)=yy(j)*dyy(j)
            fyyz(j)=yy(j)*dyz(j)
            fyzz(j)=yy(j)*dzz(j)
            fzzz(j)=zz(j)*dzz(j)

            qx(j,k+1)=qx(j,k+1)-fxxx(j)+px(j)+px(j)
            qy(j,k+1)=qy(j,k+1)-fxxy(j)
            qz(j,k+1)=qz(j,k+1)-fxxz(j)

            qx(j,k+2)=qx(j,k+2)-fxxy(j)+py(j)
            qy(j,k+2)=qy(j,k+2)-fxyy(j)+px(j)
            qz(j,k+2)=qz(j,k+2)-fxyz(j)

            qx(j,k+3)=qx(j,k+3)-fxxz(j)+pz(j)
            qy(j,k+3)=qy(j,k+3)-fxyz(j)
            qz(j,k+3)=qz(j,k+3)-fxzz(j)+px(j)

            qx(j,k+4)=qx(j,k+4)-fxyy(j)
            qy(j,k+4)=qy(j,k+4)-fyyy(j)+py(j)+py(j)
            qz(j,k+4)=qz(j,k+4)-fyyz(j)

            qx(j,k+5)=qx(j,k+5)-fxyz(j)
            qy(j,k+5)=qy(j,k+5)-fyyz(j)+pz(j)
            qz(j,k+5)=qz(j,k+5)-fyzz(j)+py(j)

            qx(j,k+6)=qx(j,k+6)-fxzz(j)
            qy(j,k+6)=qy(j,k+6)-fyzz(j)
            qz(j,k+6)=qz(j,k+6)-fzzz(j)+pz(j)+pz(j)

            fxxx(j)=fxxx(j)*twozeta
            fxxy(j)=fxxy(j)*twozeta
            fxxz(j)=fxxz(j)*twozeta
            fxyy(j)=fxyy(j)*twozeta
            fxyz(j)=fxyz(j)*twozeta
            fxzz(j)=fxzz(j)*twozeta
            fyyy(j)=fyyy(j)*twozeta
            fyyz(j)=fyyz(j)*twozeta
            fyzz(j)=fyzz(j)*twozeta
            fzzz(j)=fzzz(j)*twozeta

            gxxxx(j)=xx(j)*fxxx(j)
            gxxxy(j)=xx(j)*fxxy(j)
            gxxxz(j)=xx(j)*fxxz(j)
            gxxyy(j)=xx(j)*fxyy(j)
            gxxyz(j)=xx(j)*fxyz(j)
            gxxzz(j)=xx(j)*fxzz(j)
            gxyyy(j)=xx(j)*fyyy(j)
            gxyyz(j)=xx(j)*fyyz(j)
            gxyzz(j)=xx(j)*fyzz(j)
            gxzzz(j)=xx(j)*fzzz(j)
            gyyyy(j)=yy(j)*fyyy(j)
            gyyyz(j)=yy(j)*fyyz(j)
            gyyzz(j)=yy(j)*fyzz(j)
            gyzzz(j)=yy(j)*fzzz(j)
            gzzzz(j)=zz(j)*fzzz(j)

            q1(j,k+1)=q1(j,k+1)+gxxxx(j)+s(j)+s(j)-dxx(j)*five
            q2(j,k+1)=q2(j,k+1)+gxxxy(j)          -dxy(j)-dxy(j)
            q3(j,k+1)=q3(j,k+1)+gxxxz(j)          -dxz(j)-dxz(j)
            q4(j,k+1)=q4(j,k+1)+gxxyy(j)          -dxx(j)
            q5(j,k+1)=q5(j,k+1)+gxxyz(j)
            q6(j,k+1)=q6(j,k+1)+gxxzz(j)          -dxx(j)

            q1(j,k+2)=q1(j,k+2)+gxxxy(j)          -dxy(j)-dxy(j)-dxy(j)
            q2(j,k+2)=q2(j,k+2)+gxxyy(j)+s(j)     -dyy(j)-dxx(j)
            q3(j,k+2)=q3(j,k+2)+gxxyz(j)          -dyz(j)
            q4(j,k+2)=q4(j,k+2)+gxyyy(j)          -dxy(j)-dxy(j)-dxy(j)
            q5(j,k+2)=q5(j,k+2)+gxyyz(j)          -dxz(j)
            q6(j,k+2)=q6(j,k+2)+gxyzz(j)          -dxy(j)

            q1(j,k+3)=q1(j,k+3)+gxxxz(j)          -dxz(j)-dxz(j)-dxz(j)
            q2(j,k+3)=q2(j,k+3)+gxxyz(j)          -dyz(j)
            q3(j,k+3)=q3(j,k+3)+gxxzz(j)+s(j)     -dzz(j)-dxx(j)
            q4(j,k+3)=q4(j,k+3)+gxyyz(j)          -dxz(j)
            q5(j,k+3)=q5(j,k+3)+gxyzz(j)          -dxy(j)
            q6(j,k+3)=q6(j,k+3)+gxzzz(j)          -dxz(j)-dxz(j)-dxz(j)

            q1(j,k+4)=q1(j,k+4)+gxxyy(j)          -dyy(j)
            q2(j,k+4)=q2(j,k+4)+gxyyy(j)          -dxy(j)-dxy(j)
            q3(j,k+4)=q3(j,k+4)+gxyyz(j)
            q4(j,k+4)=q4(j,k+4)+gyyyy(j)+s(j)+s(j)-dyy(j)*five
            q5(j,k+4)=q5(j,k+4)+gyyyz(j)          -dyz(j)-dyz(j)
            q6(j,k+4)=q6(j,k+4)+gyyzz(j)          -dyy(j)

            q1(j,k+5)=q1(j,k+5)+gxxyz(j)          -dyz(j)
            q2(j,k+5)=q2(j,k+5)+gxyyz(j)          -dxz(j)
            q3(j,k+5)=q3(j,k+5)+gxyzz(j)          -dxy(j)
            q4(j,k+5)=q4(j,k+5)+gyyyz(j)          -dyz(j)-dyz(j)-dyz(j)
            q5(j,k+5)=q5(j,k+5)+gyyzz(j)+s(j)     -dzz(j)-dyy(j)
            q6(j,k+5)=q6(j,k+5)+gyzzz(j)          -dyz(j)-dyz(j)-dyz(j)

            q1(j,k+6)=q1(j,k+6)+gxxzz(j)          -dzz(j)
            q2(j,k+6)=q2(j,k+6)+gxyzz(j)
            q3(j,k+6)=q3(j,k+6)+gxzzz(j)          -dxz(j)-dxz(j)
            q4(j,k+6)=q4(j,k+6)+gyyzz(j)          -dzz(j)
            q5(j,k+6)=q5(j,k+6)+gyzzz(j)          -dyz(j)-dyz(j)
            q6(j,k+6)=q6(j,k+6)+gzzzz(j)+s(j)+s(j)-dzz(j)*five

                               endif
 1022 continue

 1020 continue

      do 1023 j=1,npts

      qf(j,k+1)=qf(j,k+1)*rsqrt3
      qf(j,k+4)=qf(j,k+4)*rsqrt3
      qf(j,k+6)=qf(j,k+6)*rsqrt3

      qx(j,k+1)=qx(j,k+1)*rsqrt3
      qy(j,k+1)=qy(j,k+1)*rsqrt3
      qz(j,k+1)=qz(j,k+1)*rsqrt3

      qx(j,k+4)=qx(j,k+4)*rsqrt3
      qy(j,k+4)=qy(j,k+4)*rsqrt3
      qz(j,k+4)=qz(j,k+4)*rsqrt3

      qx(j,k+6)=qx(j,k+6)*rsqrt3
      qy(j,k+6)=qy(j,k+6)*rsqrt3
      qz(j,k+6)=qz(j,k+6)*rsqrt3

      q1(j,k+1)=q1(j,k+1)*rsqrt3
      q2(j,k+1)=q2(j,k+1)*rsqrt3
      q3(j,k+1)=q3(j,k+1)*rsqrt3
      q4(j,k+1)=q4(j,k+1)*rsqrt3
      q5(j,k+1)=q5(j,k+1)*rsqrt3
      q6(j,k+1)=q6(j,k+1)*rsqrt3

      q1(j,k+4)=q1(j,k+4)*rsqrt3
      q2(j,k+4)=q2(j,k+4)*rsqrt3
      q3(j,k+4)=q3(j,k+4)*rsqrt3
      q4(j,k+4)=q4(j,k+4)*rsqrt3
      q5(j,k+4)=q5(j,k+4)*rsqrt3
      q6(j,k+4)=q6(j,k+4)*rsqrt3

      q1(j,k+6)=q1(j,k+6)*rsqrt3
      q2(j,k+6)=q2(j,k+6)*rsqrt3
      q3(j,k+6)=q3(j,k+6)*rsqrt3
      q4(j,k+6)=q4(j,k+6)*rsqrt3
      q5(j,k+6)=q5(j,k+6)*rsqrt3
      q6(j,k+6)=q6(j,k+6)*rsqrt3

 1023 continue

      extent3(k+1)=sumsum(npts,qf(1,k+1))/real64(npts)
      extent3(k+2)=sumsum(npts,qf(1,k+2))/real64(npts)
      extent3(k+3)=sumsum(npts,qf(1,k+3))/real64(npts)
      extent3(k+4)=sumsum(npts,qf(1,k+4))/real64(npts)
      extent3(k+5)=sumsum(npts,qf(1,k+5))/real64(npts)
      extent3(k+6)=sumsum(npts,qf(1,k+6))/real64(npts)

 1016 k=k+6

      return
      end
