      subroutine elf_dif(elfd,delf,xp,yp,zp)
      use orbitals_module
      use confs_module
      use eval_module
      implicit none
      double precision rha,rhb,gradfa,gradfb,gradf,dgradf(3)
      double precision drha(3),drhb(3)
      double precision dgradfa(3),dgradfb(3)
      double precision dratio(3),p163
      double precision ratio,x,y,z,r2
      double precision x0,x1,x2,y0,y1,y2,z0,z1,z2
      double precision alfa2,alfax,alfax2,alfay,alfay2,alfaz,alfaz2
      double precision zp,yp,xp,ex,fermi,tol
      double precision epsex
      double precision rho,elfd,drho(3),delf(3)
      integer i, j, k, l, mu, level,imu
      integer ij,kl,imol,jmol,lmol,kmol
      common/accuracy/level
      data fermi/0.121300565d0/
      p163=-16.0d0/3.0d0
      tol=2.871d-5
      epsex=0.1d0**(7-level)
      if(ifci) then
         nmola=nmol
      endif
      do i=1,nmola
         indij(i)=(i-1)*(nmola+nmola-i)/2
      enddo
      imu=0
      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
            if(i.eq.0) then
              x2=alfa2*(1.0d0+alfax2)
            else if (i.eq.1) then
              x2=alfax*(3.d0+alfax2)
            else if (i.eq.2) then
              x2=2.0d0+alfax2*(5.0d0+alfax2)
            else if (i.eq.3) then
              x2=x*(6.0d0+alfax2*(7.d0+alfax2))
            else if (i.gt.3) then
              x2=x**(i-2)*(dfloat(i)*dfloat(i-1)+alfax2*
     .          (dfloat(i+i+1)+alfax2))
            endif
            if(j.eq.0) then
               y2=alfa2*(1.0d0+alfay2)
            else if (j.eq.1) then
               y2=alfay*(3.d0+alfay2)
            else if (j.eq.2) then
               y2=2.0d0+alfay2*(5.0d0+alfay2)
            else if (j.eq.3) then
               y2=y*(6.0d0+alfay2*(7.d0+alfay2))
            else if (j.gt.3) then
               y2=y**(j-2)*(dfloat(j)*dfloat(j-1)+alfay2*
     .            (dfloat(j+j+1)+alfay2))
            endif
            if(k.eq.0) then
               z2=alfa2*(1.0d0+alfaz2)
            else if (k.eq.1) then
               z2=alfaz*(3.d0+alfaz2)
            else if (k.eq.2) then
               z2=2.0d0+alfaz2*(5.0d0+alfaz2)
            else if (k.eq.3) then
               z2=z*(6.0d0+alfaz2*(7.d0+alfaz2))
            else if (k.gt.3) then
               z2=z**(k-2)*(dfloat(k)*dfloat(k-1)+alfaz2*
     .            (dfloat(k+k+1)+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
            fxx(imu)=x2*y0*z0*ex
            fxy(imu)=x1*y1*z0*ex
            fxz(imu)=x1*y0*z1*ex
            fyy(imu)=x0*y2*z0*ex
            fyz(imu)=x0*y1*z1*ex
            fzz(imu)=x0*y0*z2*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
         d2phi(imol,1)=0.0d0
         d2phi(imol,2)=0.0d0
         d2phi(imol,3)=0.0d0
         d2phi(imol,4)=0.0d0
         d2phi(imol,5)=0.0d0
         d2phi(imol,6)=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)
            d2phi(imol,1)=d2phi(imol,1)+c(imol,mu)*fxx(i)
            d2phi(imol,2)=d2phi(imol,2)+c(imol,mu)*fxy(i)
            d2phi(imol,3)=d2phi(imol,3)+c(imol,mu)*fxz(i)
            d2phi(imol,4)=d2phi(imol,4)+c(imol,mu)*fyy(i)
            d2phi(imol,5)=d2phi(imol,5)+c(imol,mu)*fyz(i)
            d2phi(imol,6)=d2phi(imol,6)+c(imol,mu)*fzz(i)
         enddo
      enddo
      do imol=1,nmola
         do jmol=imol,nmola
            lmol=indij(imol)+jmol
            rhorb(lmol)=phi(imol)*phi(jmol)
            dotorb(lmol,1)=phi(imol)*dphi(jmol,1)- 
     .                  phi(jmol)*dphi(imol,1)
            dotorb(lmol,2)=phi(imol)*dphi(jmol,2)- 
     .                  phi(jmol)*dphi(imol,2)
            dotorb(lmol,3)=phi(imol)*dphi(jmol,3)- 
     .                  phi(jmol)*dphi(imol,3)
            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)
            ddotorb(lmol,1)=phi(imol)*d2phi(jmol,1)-
     .                      phi(jmol)*d2phi(imol,1)
            ddotorb(lmol,2)=phi(imol)*d2phi(jmol,2)+
     .                      dphi(imol,2)*dphi(jmol,1)-
     .                      phi(jmol)*d2phi(imol,2)-
     .                      dphi(imol,1)*dphi(jmol,2)
            ddotorb(lmol,3)=phi(imol)*d2phi(jmol,3)+
     .                      dphi(imol,3)*dphi(jmol,1)-
     .                      phi(jmol)*d2phi(imol,3)-
     .                      dphi(imol,1)*dphi(jmol,3)
            ddotorb(lmol,4)=phi(imol)*d2phi(jmol,2)+
     .                      dphi(imol,1)*dphi(jmol,2)-
     .                      phi(jmol)*d2phi(imol,2)-
     .                      dphi(imol,2)*dphi(jmol,1)
            ddotorb(lmol,5)=phi(imol)*d2phi(jmol,4)- 
     .                      phi(jmol)*d2phi(imol,4)
            ddotorb(lmol,6)=phi(imol)*d2phi(jmol,5)+
     .                      dphi(imol,3)*dphi(jmol,2)-
     .                      phi(jmol)*d2phi(imol,5)-
     .                      dphi(imol,2)*dphi(jmol,3)
            ddotorb(lmol,7)=phi(imol)*d2phi(jmol,3)+
     .                      dphi(imol,1)*dphi(jmol,3)-
     .                      phi(jmol)*d2phi(imol,3)-
     .                      dphi(imol,3)*dphi(jmol,1)
            ddotorb(lmol,8)=phi(imol)*d2phi(jmol,5)+
     .                      dphi(imol,2)*dphi(jmol,3)-
     .                      phi(jmol)*d2phi(imol,5)-
     .                      dphi(imol,3)*dphi(jmol,2)
            ddotorb(lmol,9)=phi(imol)*d2phi(jmol,6)- 
     .                      phi(jmol)*d2phi(imol,6)
            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
      gradfa=0.0d0
      dgradfa(1)=0.0d0
      dgradfa(2)=0.0d0
      dgradfa(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
*              compute D(r)
      do l=1,naa
         imol=iaa(l)
         jmol=jaa(l)
         kmol=kaa(l)
         lmol=laa(l)
         ij=indij(imol)+jmol
         kl=indij(kmol)+lmol
         gradfa=gradfa+paa(l)*(dotorb(ij,1)*dotorb(kl,1)+
     .          dotorb(ij,2)*dotorb(kl,2)+dotorb(ij,3)*dotorb(kl,3))
         dgradfa(1)=dgradfa(1)+paa(l)*(dotorb(ij,1)*ddotorb(kl,1)+
     .              ddotorb(ij,1)*dotorb(kl,1)+
     .              ddotorb(ij,4)*dotorb(kl,2)+
     .              dotorb(ij,2)*ddotorb(kl,4)+
     .              ddotorb(ij,7)*dotorb(kl,3)+
     .              dotorb(ij,3)*ddotorb(kl,7))
         dgradfa(2)=dgradfa(2)+paa(l)*(dotorb(ij,1)*ddotorb(kl,2)+
     .              ddotorb(ij,2)*dotorb(kl,1)+
     .              ddotorb(ij,5)*dotorb(kl,2)+
     .              dotorb(ij,2)*ddotorb(kl,5)+
     .              ddotorb(ij,8)*dotorb(kl,3)+
     .              dotorb(ij,3)*ddotorb(kl,8))
         dgradfa(3)=dgradfa(3)+paa(l)*(dotorb(ij,1)*ddotorb(kl,3)+
     .              ddotorb(ij,3)*dotorb(kl,1)+
     .              ddotorb(ij,6)*dotorb(kl,2)+
     .              dotorb(ij,2)*ddotorb(kl,6)+
     .              ddotorb(ij,9)*dotorb(kl,3)+
     .              dotorb(ij,3)*ddotorb(kl,9))
      enddo
      if(nb.eq.0) then
         rhb=rha
         gradfb=gradfa
         drhb(1)=drha(1)
         drhb(2)=drha(2)
         drhb(3)=drha(3)
         dgradfb(1)=dgradfa(1)
         dgradfb(2)=dgradfa(2)
         dgradfb(3)=dgradfa(3)
      else
      do imol=nmola+1,nmol
         do jmol=imol,nmol
            lmol=indij(imol-nmola)+jmol-nmola
            rhorb(lmol)=phi(imol)*phi(jmol)
            dotorb(lmol,1)=phi(imol)*dphi(jmol,1)-phi(jmol)*dphi(imol,1)
            dotorb(lmol,2)=phi(imol)*dphi(jmol,2)-phi(jmol)*dphi(imol,2)
            dotorb(lmol,3)=phi(imol)*dphi(jmol,3)-phi(jmol)*dphi(imol,3)
            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)
            ddotorb(lmol,1)=phi(imol)*d2phi(jmol,1)-
     .                      phi(jmol)*d2phi(imol,1)
            ddotorb(lmol,2)=phi(imol)*d2phi(jmol,2)+
     .                      dphi(imol,2)*dphi(jmol,1)-
     .                      phi(jmol)*d2phi(imol,2)-
     .                      dphi(imol,1)*dphi(jmol,2)
            ddotorb(lmol,3)=phi(imol)*d2phi(jmol,3)+
     .                      dphi(imol,3)*dphi(jmol,1)-
     .                      phi(jmol)*d2phi(imol,3)-
     .                      dphi(imol,1)*dphi(jmol,3)
            ddotorb(lmol,4)=phi(imol)*d2phi(jmol,2)+
     .                      dphi(imol,1)*dphi(jmol,2)-
     .                      phi(jmol)*d2phi(imol,2)-
     .                      dphi(imol,2)*dphi(jmol,1)
            ddotorb(lmol,5)=phi(imol)*d2phi(jmol,4)-
     .                      phi(jmol)*d2phi(imol,4)
            ddotorb(lmol,6)=phi(imol)*d2phi(jmol,5)+
     .                      dphi(imol,3)*dphi(jmol,2)-
     .                      phi(jmol)*d2phi(imol,5)-
     .                      dphi(imol,2)*dphi(jmol,3)
            ddotorb(lmol,7)=phi(imol)*d2phi(jmol,3)+
     .                      dphi(imol,1)*dphi(jmol,3)-
     .                      phi(jmol)*d2phi(imol,3)-
     .                      dphi(imol,3)*dphi(jmol,1)
            ddotorb(lmol,8)=phi(imol)*d2phi(jmol,5)+
     .                      dphi(imol,2)*dphi(jmol,3)-
     .                      phi(jmol)*d2phi(imol,5)-
     .                      dphi(imol,3)*dphi(jmol,2)
            ddotorb(lmol,9)=phi(imol)*d2phi(jmol,6)-
     .                      phi(jmol)*d2phi(imol,6)
            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
      gradfb=0.0d0
      drhb(1)=0.0d0
      drhb(2)=0.0d0
      drhb(3)=0.0d0
      dgradfb(1)=0.0d0
      dgradfb(2)=0.0d0
      dgradfb(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
      do l=1,nbb
         imol=ibb(l)
         jmol=jbb(l)
         kmol=kbb(l)
         lmol=lbb(l)
         ij=indij(imol)+jmol
         kl=indij(kmol)+lmol
         gradfb=gradfb+pbb(l)*(dotorb(ij,1)*dotorb(kl,1)+
     .          dotorb(ij,2)*dotorb(kl,2)+dotorb(ij,3)*dotorb(kl,3))
         dgradfb(1)=dgradfb(1)+pbb(l)*(dotorb(ij,1)*ddotorb(kl,1)+
     .              ddotorb(ij,1)*dotorb(kl,1)+
     .              ddotorb(ij,4)*dotorb(kl,2)+
     .              dotorb(ij,2)*ddotorb(kl,4)+
     .              ddotorb(ij,7)*dotorb(kl,3)+
     .              dotorb(ij,3)*ddotorb(kl,7))
         dgradfb(2)=dgradfb(2)+pbb(l)*(dotorb(ij,1)*ddotorb(kl,2)+
     .              ddotorb(ij,2)*dotorb(kl,1)+
     .              ddotorb(ij,5)*dotorb(kl,2)+
     .              dotorb(ij,2)*ddotorb(kl,5)+
     .              ddotorb(ij,8)*dotorb(kl,3)+
     .              dotorb(ij,3)*ddotorb(kl,8))
         dgradfb(3)=dgradfb(3)+pbb(l)*(dotorb(ij,1)*ddotorb(kl,3)+
     .              ddotorb(ij,3)*dotorb(kl,1)+
     .              ddotorb(ij,6)*dotorb(kl,2)+
     .              dotorb(ij,2)*ddotorb(kl,6)+
     .              ddotorb(ij,9)*dotorb(kl,3)+
     .              dotorb(ij,3)*ddotorb(kl,9))
      enddo
      endif
      rho=rha+rhb+1.0d-16
      gradf=(gradfa+gradfb)+tol*rho
      ratio=gradf*gradf*fermi*rho**p163
      elfd=1.d0/(1.0d0+ratio)
      drho(1)=drha(1)+drhb(1)
      drho(2)=drha(2)+drhb(2)
      drho(3)=drha(3)+drhb(3)
      dgradf(1)=dgradfa(1)+dgradfb(1)+tol*drho(1)
      dgradf(2)=dgradfa(2)+dgradfb(2)+tol*drho(2)
      dgradf(3)=dgradfa(3)+dgradfb(3)+tol*drho(3)
      dratio(1)=2.0d0*ratio*dgradf(1)/gradf+p163*ratio*drho(1)/rho
      dratio(2)=2.0d0*ratio*dgradf(2)/gradf+p163*ratio*drho(2)/rho
      dratio(3)=2.0d0*ratio*dgradf(3)/gradf+p163*ratio*drho(3)/rho
      delf(1)=-elfd*elfd*dratio(1)
      delf(2)=-elfd*elfd*dratio(2)
      delf(3)=-elfd*elfd*dratio(3)
      return
      end
