      subroutine munu_contrib(mu,nu)
      use orbitals_module
      use eval_module
      use grid_module
      implicit none
      double precision sqrpi, int_gauss(0:8), sum
      double precision x, y, z
      double precision 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
      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
      double precision dx,dy,dz
      integer ix,iy,iz,indice
      integer s, t, la, lb, lmax, mu, nu
      integer i
      common/grid/ dx,dy,dz
      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)
      xa=rc(1,mu)
      ya=rc(2,mu)
      za=rc(3,mu)
      xb=rc(1,nu)
      yb=rc(2,nu)
      zb=rc(3,nu)
      alfap=alfa(mu)+alfa(nu)
      sqralfap=dsqrt(alfap)
*
*        contribution 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=xmin-dx
      do ix=1,np(1)
         x=x+dx
         x_min=x-0.5d0*dx-xp
         x_max=x_min+dx
         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
         gx(ix)=sum*dexp(arg)
      enddo
      gx(1)=0.5d0*gx(1)
      gx(np(1))=0.5d0*gx(np(1))
*
*        contribution 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=ymin-dy
      do iy=1,np(2)
         y=y+dy
         y_min=y-0.5d0*dy-yp
         y_max=y_min+dy
         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
         gy(iy)=sum*dexp(arg)
      enddo
      gy(1)=0.5d0*gy(1)
      gy(np(2))=0.5d0*gy(np(2))
*
*        contribution 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=zmin-dz
      do iz=1,np(3)
         z=z+dz
         z_min=z-0.5d0*dz-zp
         z_max=z_min+dz
         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
         gz(iz)=sum*dexp(arg)
      enddo
      gz(1)=0.5d0*gz(1)
      gz(np(3))=0.5d0*gz(np(3))
      return
      end
