      subroutine rho_dif(rho,drho,xp,yp,zp)
      use grid_module
      use confs_module
      use orbitals_module
      use eval_module
      implicit none
      double precision x,y,z
      double precision x0,x1,y0,y1,z0,z1,r2
      double precision alfa2,alfax,alfax2,alfay,alfay2,alfaz,alfaz2
      double precision zp,yp,xp,ex,tol
      double precision epsex
      double precision rho,drho(3)
      double precision rha,rhb,drha(3),drhb(3)
      integer i, j, k, l, mu, level,imu
      integer imol,jmol,lmol
      common/threshold/ tol
      common/accuracy/level
      tol=2.871d-5
      epsex=0.1d0**(7-level)
      imu=0
      if(ifci) then
         nmola=nmol
      endif
         
      do i=1,nmola
         indij(i)=(i-1)*(nmola+nmola-i)/2
      enddo
      do mu=1, ngto
         x=xp-rc(1,mu)
         y=yp-rc(2,mu)
         z=zp-rc(3,mu)
         r2=x*x+y*y+z*z
         ex=dexp(-alfa(mu)*r2)
         if(ex.gt.epsex) then
            imu=imu+1
            extest(imu)=mu
            alfa2=-alfa(mu)*2.0d0
            alfax=alfa2*x
            alfax2=alfax*x
            alfay=alfa2*y
            alfay2=alfay*y
            alfaz=alfa2*z
            alfaz2=alfaz*z
            i=ir(1,mu)
            j=ir(2,mu)
            k=ir(3,mu)
            if(i.eq.0) then
               x0=1.0d0
               x1=alfax
            else if (i.eq.1) then
               x0=x
               x1=1.0d0+alfax2
            else if (i.eq.2) then
               x0=x*x
               x1=x*(2.d0+alfax2)
            else if (i.eq.3) then
               x0=x**3
               x1=x*x*(3.0d0+alfax2)
            else if (i.gt.3) then
               x0=x**i
               x1=x**(i-1)*(dfloat(i)+alfax2)
            endif
            if(j.eq.0) then
               y0=1.0d0
               y1=alfay
            else if (j.eq.1) then
               y0=y
               y1=1.0d0+alfay2
            else if (j.eq.2) then
               y0=y*y
               y1=y*(2.d0+alfay2)
            else if (j.eq.3) then
               y0=y**3
               y1=y*y*(3.0d0+alfay2)
            else if (j.gt.3) then
               y0=y**j
               y1=y**(j-1)*(dfloat(j)+alfay2)
            endif
            if(k.eq.0) then
               z0=1.0d0
               z1=alfaz
            else if (k.eq.1) then
               z0=z
               z1=1.0d0+alfaz2
            else if (k.eq.2) then
               z0=z*z
               z1=z*(2.d0+alfaz2)
            else if (k.eq.3) then
               z0=z**3
               z1=z*z*(3.0d0+alfaz2)
            else if (k.gt.3) then
               z0=z**k
               z1=z**(k-1)*(dfloat(k)+alfaz2)
            endif
            g(imu)=x0*y0*z0*ex
            gx(imu)=x1*y0*z0*ex
            gy(imu)=x0*y1*z0*ex
            gz(imu)=x0*y0*z1*ex
         endif
      enddo
      do imol=1,nmol
         phi(imol)=0.0d0
         dphi(imol,1)=0.0d0
         dphi(imol,2)=0.0d0
         dphi(imol,3)=0.0d0
         do i=1,imu
            mu=extest(i)
            phi(imol)=phi(imol)+c(imol,mu)*g(i)
            dphi(imol,1)=dphi(imol,1)+c(imol,mu)*gx(i)
            dphi(imol,2)=dphi(imol,2)+c(imol,mu)*gy(i)
            dphi(imol,3)=dphi(imol,3)+c(imol,mu)*gz(i)
         enddo
      enddo
      do imol=1,nmola
         do jmol=imol,nmola
            lmol=indij(imol)+jmol
            rhorb(lmol)=phi(imol)*phi(jmol)
            drhorb(lmol,1)=phi(imol)*dphi(jmol,1)+
     .                     phi(jmol)*dphi(imol,1)
            drhorb(lmol,2)=phi(imol)*dphi(jmol,2)+
     .                     phi(jmol)*dphi(imol,2)
            drhorb(lmol,3)=phi(imol)*dphi(jmol,3)+
     .                     phi(jmol)*dphi(imol,3)
            if(imol.ne.jmol) then
                rhorb(lmol)=rhorb(lmol)+rhorb(lmol)
                drhorb(lmol,1)=drhorb(lmol,1)+drhorb(lmol,1)
                drhorb(lmol,2)=drhorb(lmol,2)+drhorb(lmol,2)
                drhorb(lmol,3)=drhorb(lmol,3)+drhorb(lmol,3)
            endif
         enddo
      enddo
      rha=0.0d0
      drha(1)=0.0d0
      drha(2)=0.0d0
      drha(3)=0.0d0
      do l=1,na
         lmol=ia(l)
         rha=rha+pa(l)*rhorb(lmol)
         drha(1)=drha(1)+pa(l)*drhorb(lmol,1)
         drha(2)=drha(2)+pa(l)*drhorb(lmol,2)
         drha(3)=drha(3)+pa(l)*drhorb(lmol,3)
      enddo
      if(nb.eq.0) then
         rhb=rha
         drhb(1)=drha(1)
         drhb(2)=drha(2)
         drhb(3)=drha(3)
      else
         do imol=nmola+1,nmol
            do jmol=imol,nmol
               lmol=indij(imol-nmola)+jmol-nmola
               rhorb(lmol)=phi(imol)*phi(jmol)
               drhorb(lmol,1)=phi(imol)*dphi(jmol,1)+
     .                        phi(jmol)*dphi(imol,1)
               drhorb(lmol,2)=phi(imol)*dphi(jmol,2)+
     .                        phi(jmol)*dphi(imol,2)
               drhorb(lmol,3)=phi(imol)*dphi(jmol,3)+
     .                        phi(jmol)*dphi(imol,3)
               if(imol.ne.jmol) then
                   rhorb(lmol)=rhorb(lmol)+rhorb(lmol)
                   drhorb(lmol,1)=drhorb(lmol,1)+drhorb(lmol,1)
                   drhorb(lmol,2)=drhorb(lmol,2)+drhorb(lmol,2)
                   drhorb(lmol,3)=drhorb(lmol,3)+drhorb(lmol,3)
               endif
            enddo
         enddo
         rhb=0.0d0
         drhb(1)=0.0d0
         drhb(2)=0.0d0
         drhb(3)=0.0d0
         do l=1,nb
            lmol=ib(l)
            rhb=rhb+pb(l)*rhorb(lmol)
            drhb(1)=drhb(1)+pb(l)*drhorb(lmol,1)
            drhb(2)=drhb(2)+pb(l)*drhorb(lmol,2)
            drhb(3)=drhb(3)+pb(l)*drhorb(lmol,3)
         enddo
      endif
      rho=rha+rhb
      drho(1)=drha(1)+drhb(1)
      drho(2)=drha(2)+drhb(2)
      drho(3)=drha(3)+drhb(3)
      return
      end
