      subroutine box_dens(rho_tol)
      use orbitals_module
      use grid_module
      use popul_module
      implicit none
      double precision sqrpi, int_gauss(0:8), sum
      double precision x, y, z, dv, estimate, r2, xx, yy, zz
      double precision sum_xyz, arg_erf, sqralfap, arg, pi
      double precision xapower, xbpower, xa, xb, xp, za, zb, ya
      double precision yapower, ybpower, zapower, zbpower
      double precision yb, alfap, xx1, yy1, zz1
      double precision binom(0:3,0:3), x_min, x_max, y_min, y_max
      double precision erf0(0:10000),erf1(0:10000),darg
      double precision z_min, z_max, yp, zp,rho_tol
      integer s, t, la, lb, lmax, mu, nu, munu
      integer i, indice
      common /erftab/ erf0,erf1
      data binom
     .   / 1.0d0,  0.0d0,  0.0d0,  0.0d0,
     .     1.0d0,  1.0d0,  0.0d0,  0.0d0,
     .     1.0d0,  2.0d0,  1.0d0,  0.0d0,
     .     1.0d0,  3.0d0,  3.0d0,  1.0d0/
      pi=dacos(-1.0d0)
      sqrpi=dsqrt(pi)
      x=0.5d0*(xmax+xmin)
      y=0.5d0*(ymax+ymin)
      z=0.5d0*(zmax+zmin)
      dv=(xmax-xmin)*(ymax-ymin)*(zmax-zmin)
      munu=0
      do mu=1,ngto
         do nu=mu,ngto
         munu=munu+1
         smunu(munu)=0.0d0
         xa=rc(1,mu)
         ya=rc(2,mu)
         za=rc(3,mu)
         xb=rc(1,nu)
         yb=rc(2,nu)
         zb=rc(3,nu)
*
*        mu nu contribution to total density
*
         xx=x-xa
         yy=y-ya
         zz=z-za
         r2= xx*xx+yy*yy+zz*zz
         xx1=1.0d0
         yy1=1.0d0
         zz1=1.0d0
         if(ir(1,mu).gt.0) then
            xx1=xx**ir(1,mu)
         endif
         if(ir(2,mu).gt.0) then
            yy1=yy**ir(2,mu)
         endif
         if(ir(3,mu).gt.0) then
            zz1=zz**ir(3,mu)
         endif
         estimate= xx1*yy1*zz1*dexp(-alfa(mu)*r2)
         xx=x-xb
         yy=y-yb
         zz=z-zb
         r2= xx*xx+yy*yy+zz*zz
         xx1=1.0d0
         yy1=1.0d0
         zz1=1.0d0
         if(ir(1,nu).gt.0) then
            xx1=xx**ir(1,nu)
         endif
         if(ir(2,nu).gt.0) then
            yy1=yy**ir(2,nu)
         endif
         if(ir(3,nu).gt.0) then
            zz1=zz**ir(3,nu)
         endif
         estimate= estimate*xx1*yy1*zz1*dexp(-alfa(nu)*r2)
*
         if(dabs(estimate).ge.rho_tol) then
         alfap=alfa(mu)+alfa(nu)
         sqralfap=dsqrt(alfap)
*
*        integration along x direction
*
         la=ir(1,mu)
         lb=ir(1,nu)
         xp=(alfa(mu)*xa+alfa(nu)*xb)/alfap
         xa=xp-xa
         xb=xp-xb
         lmax=la+lb
         x_min=xmin-xp
         x_max=xmax-xp
         arg_erf=sqralfap*dabs(x_min)
         indice=nint(2.0d3*arg_erf)
         if(indice.le.10000) then
            darg=arg_erf-5.0d-4*dfloat(indice)
            int_gauss(0)=erf0(indice)+darg*erf1(indice)
         else
            int_gauss(0)=1.0d0
         endif
         if(x_min.lt.0.0d0) then
            int_gauss(0)=-int_gauss(0)
         endif
         arg_erf=sqralfap*dabs(x_max)
         if(x_max.lt.0.0d0) then
            indice=nint(2.0d3*arg_erf)
            if(indice.le.10000) then
               darg=arg_erf-5.0d-4*dfloat(indice)
               int_gauss(0)=-erf0(indice)-erf1(indice)*darg-int_gauss(0)
            else
               int_gauss(0)=-1.0d0-int_gauss(0)
            endif
         else
            indice=nint(2.0d3*arg_erf)
            if(indice.le.10000) then
               darg=arg_erf-5.0d-4*dfloat(indice)
               int_gauss(0)=erf0(indice)+erf1(indice)*darg-int_gauss(0)
            else
               int_gauss(0)=1.0d0-int_gauss(0)
            endif
         endif
         int_gauss(1)=0.5d0*(dexp(-alfap*x_min*x_min)-
     .                dexp(-alfap*x_max*x_max))/alfap
         int_gauss(0)=0.5d0*sqrpi*int_gauss(0)/sqralfap
         do i=2,lmax
            int_gauss(i)=0.5d0*dfloat(i-1)*int_gauss(i-2)/alfap+
     .                   0.5d0*(x_min**(i-1)*dexp(-alfap*x_min*x_min)-
     .                   x_max**(i-1)*dexp(-alfap*x_max*x_max))/alfap
         enddo
         xapower=1.0d0
         sum=0.0d0
         do s=0,la
            xbpower=1.0d0
            do t=0,lb
               sum=sum+binom(s,la)*binom(t,lb)*xapower*xbpower*
     .             int_gauss(lmax-s-t)
               xbpower=xbpower*xb
            enddo
         xapower=xapower*xa
         enddo
         arg=-alfa(mu)*alfa(nu)*(xa-xb)*(xa-xb)/alfap
         sum_xyz=sum*dexp(arg)
*
*        integration along y direction
*
         la=ir(2,mu)
         lb=ir(2,nu)
         yp=(alfa(mu)*ya+alfa(nu)*yb)/alfap
         ya=yp-ya
         yb=yp-yb
         lmax=la+lb
         y_min=ymin-yp
         y_max=ymax-yp
         arg_erf=sqralfap*dabs(y_min)
         indice=nint(2.0d3*arg_erf)
         if(indice.le.10000) then
            darg=arg_erf-5.0d-4*dfloat(indice)
            int_gauss(0)=erf0(indice)+darg*erf1(indice)
         else
            int_gauss(0)=1.0d0
         endif
         if(y_min.lt.0.0d0) then
            int_gauss(0)=-int_gauss(0)
         endif
         arg_erf=sqralfap*dabs(y_max)
         if(y_max.lt.0.0d0) then
            indice=nint(2.0d3*arg_erf)
            if(indice.le.10000) then
               darg=arg_erf-5.0d-4*dfloat(indice)
               int_gauss(0)=-erf0(indice)-erf1(indice)*darg-int_gauss(0)
            else
               int_gauss(0)=-1.0d0-int_gauss(0)
            endif
         else
           indice=nint(2.0d3*arg_erf)
            if(indice.le.10000) then
               darg=arg_erf-5.0d-4*dfloat(indice)
               int_gauss(0)=erf0(indice)+erf1(indice)*darg-int_gauss(0)
            else
               int_gauss(0)=1.0d0-int_gauss(0)
            endif
         endif
         int_gauss(0)=0.5d0*sqrpi*int_gauss(0)/sqralfap
         int_gauss(1)=0.5d0*(dexp(-alfap*y_min*y_min)-
     .                dexp(-alfap*y_max*y_max))/alfap
         do i=2,lmax
            int_gauss(i)=0.5d0*dfloat(i-1)*int_gauss(i-2)/alfap+
     .                   0.5d0*(y_min**(i-1)*dexp(-alfap*y_min*y_min)-
     .                   y_max**(i-1)*dexp(-alfap*y_max*y_max))/alfap
         enddo
         yapower=1.0d0
         sum=0.0d0
         do s=0,la
            ybpower=1.0d0
            do t=0,lb
               sum=sum+binom(s,la)*binom(t,lb)*yapower*ybpower*
     .             int_gauss(lmax-s-t)
               ybpower=ybpower*yb
            enddo
         yapower=yapower*ya
         enddo
         arg=-alfa(mu)*alfa(nu)*(ya-yb)*(ya-yb)/alfap
         sum_xyz=sum_xyz*sum*dexp(arg)
*
*        integration along z direction
*
         la=ir(3,mu)
         lb=ir(3,nu)
         zp=(alfa(mu)*za+alfa(nu)*zb)/alfap
         za=zp-za
         zb=zp-zb
         lmax=la+lb
         z_min=zmin-zp
         z_max=zmax-zp
         arg_erf=sqralfap*dabs(z_min)
         indice=nint(2.0d3*arg_erf)
         if(indice.le.10000) then
            darg=arg_erf-5.0d-4*dfloat(indice)
            int_gauss(0)=erf0(indice)+darg*erf1(indice)
         else
            int_gauss(0)=1.0d0
         endif
         if(z_min.lt.0.0d0) then
            int_gauss(0)=-int_gauss(0)
         endif
         arg_erf=sqralfap*dabs(z_max)
         if(z_max.lt.0.0d0) then
            indice=nint(2.0d3*arg_erf)
            if(indice.le.10000) then
               darg=arg_erf-5.0d-4*dfloat(indice)
               int_gauss(0)=-erf0(indice)-erf1(indice)*darg-int_gauss(0)
            else
               int_gauss(0)=-1.0d0-int_gauss(0)
            endif
         else
           indice=nint(2.0d3*arg_erf)
            if(indice.le.10000) then
               darg=arg_erf-5.0d-4*dfloat(indice)
               int_gauss(0)=erf0(indice)+erf1(indice)*darg-int_gauss(0)
            else
               int_gauss(0)=1.0d0-int_gauss(0)
            endif
         endif
         int_gauss(0)=0.5d0*sqrpi*int_gauss(0)/sqralfap
         int_gauss(1)=0.5d0*(dexp(-alfap*z_min*z_min)-
     .                dexp(-alfap*z_max*z_max))/alfap
         do i=2,lmax
            int_gauss(i)=0.5d0*dfloat(i-1)*int_gauss(i-2)/alfap+
     .                   0.5d0*(z_min**(i-1)*dexp(-alfap*z_min*z_min)-
     .                   z_max**(i-1)*dexp(-alfap*z_max*z_max))/alfap
         enddo
         zapower=1.0d0
         sum=0.0d0
         do s=0,la
            zbpower=1.0d0
            do t=0,lb
               sum=sum+binom(s,la)*binom(t,lb)*zapower*zbpower*
     .             int_gauss(lmax-s-t)
               zbpower=zbpower*zb
            enddo
         zapower=zapower*za
         enddo
         arg=-alfa(mu)*alfa(nu)*(za-zb)*(za-zb)/alfap
         sum_xyz=sum_xyz*sum*dexp(arg)
         else
         sum_xyz=estimate*dv
         endif
         smunu(munu)=smunu(munu)+sum_xyz
      enddo
      enddo
      return
      end
