      subroutine compg1(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,
     &                  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)

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 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(*)

      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
 1003 qz(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)

                               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
      qf(j,k+2)=zero
      qx(j,k+2)=zero
      qy(j,k+2)=zero
      qz(j,k+2)=zero
      qf(j,k+3)=zero
      qx(j,k+3)=zero
      qy(j,k+3)=zero
 1010 qz(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)

                               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
      qf(j,k+2)=zero
      qx(j,k+2)=zero
      qy(j,k+2)=zero
      qz(j,k+2)=zero
      qf(j,k+3)=zero
      qx(j,k+3)=zero
      qy(j,k+3)=zero
      qz(j,k+3)=zero
      qf(j,k+4)=zero
      qx(j,k+4)=zero
      qy(j,k+4)=zero
      qz(j,k+4)=zero
      qf(j,k+5)=zero
      qx(j,k+5)=zero
      qy(j,k+5)=zero
      qz(j,k+5)=zero
      qf(j,k+6)=zero
      qx(j,k+6)=zero
      qy(j,k+6)=zero
 1017 qz(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)

                               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

 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
