      subroutine elfcalc(thresh)
      use orbitals_module
      use grid_module
      use confs_module
      implicit none
      double precision rh,rha,rhb,gradfa,gradfb
      double precision gradf,fermi,tol,elfd,thresh
      double precision ratio
      double precision lapl,lapla,laplb
      double precision, dimension(:),allocatable :: d, dx, dy,dz
      double precision, dimension(:),allocatable :: dx2, dy2,dz2
      double precision, dimension(:),allocatable :: phi,rhorb,laporb
      double precision, dimension(:,:),allocatable :: dphi,d2phi,dotorb
      integer i, j, k, l,  mu, nu, munu, ibuf
      integer iaux, jaux, kaux, jbuf
      integer imol,jmol,kmol,lmol,ij,kl
      data fermi/0.121300565d0/
      allocate(phi(nmol),dphi(nmol,3),d2phi(nmol,3),rhorb(nmol*
     .         (nmol+1)/2))
      allocate(laporb(nmol*(nmol+1)/2),dotorb(nmol*(nmol+1)/2,3))
      allocate(d(ngto),dx(ngto),dy(ngto),dz(ngto),dx2(ngto),dy2(ngto),
     .         dz2(ngto))
!
      ibuf=np(1)*np(2)*np(3)
      do i=1,ibuf
         buffer_elf(i)=0.0e0
         buffer_rho(i)=0.0e0
      enddo
      tol=2.871d-5
!     grid point selection
      do i=1,np(1),4
         iaux=i
         if(i.gt.1) then
            iaux=i-1
         endif
         do j=1,np(2),4
            jaux=j
            if(j.gt.1) then
               jaux=j-1
            endif
            do k=1,np(3),4
               kaux=k
               if(k.gt.1) then
                  kaux=k-1
               endif
               rh=0.0d0
               do mu=1, ngto
                  d(mu)=fx(i,mu)*fy(j,mu)*fz(k,mu)
               enddo
               munu=0
               do mu=1,ngto
                  if(dabs(d(mu)).gt.1.d-5) then
                     do nu=mu,ngto
                       munu=munu+1 
                       if(dabs(d(nu)).gt.1.d-5) then
                          rh=rh+p(munu)*d(mu)*d(nu)
                       endif
                     enddo
                  else
                     munu=munu+ngto+1-mu
                  endif
               enddo
               if(rh.gt.thresh) then
                  ibuf=k+np(3)*(j-1+np(2)*(i-1))
                  buffer_rho(ibuf)=sngl(rh)
                  ibuf=kaux+np(3)*(j-1+np(2)*(i-1))
                  buffer_rho(ibuf)=sngl(rh)
                  ibuf=k+np(3)*(jaux-1+np(2)*(i-1))
                  buffer_rho(ibuf)=sngl(rh)
                  ibuf=kaux+np(3)*(jaux-1+np(2)*(i-1))
                  buffer_rho(ibuf)=sngl(rh)
                  ibuf=k+np(3)*(j-1+np(2)*(iaux-1))
                  buffer_rho(ibuf)=sngl(rh)
                  ibuf=kaux+np(3)*(j-1+np(2)*(iaux-1))
                  buffer_rho(ibuf)=sngl(rh)
                  ibuf=k+np(3)*(jaux-1+np(2)*(iaux-1))
                  buffer_rho(ibuf)=sngl(rh)
                  ibuf=kaux+np(3)*(jaux-1+np(2)*(iaux-1))
                  buffer_rho(ibuf)=sngl(rh)
               endif
            enddo
         enddo
      enddo
!     end  selection 
      do i=1,np(1)
         iaux=4*int(i/4)
         if(iaux.eq.0) then
            iaux=1
         endif
         do j=1,np(2)
            jaux=4*int(j/4)
            if(jaux.eq.0) then
               jaux=1
            endif
            do k=1,np(3)
               kaux=1+4*int(k/4)
               if(kaux.eq.0) then
                  kaux=1
               endif
               ibuf=k+np(3)*(j-1+np(2)*(i-1))
               jbuf=kaux+np(3)*(jaux-1+np(2)*(iaux-1))
               if(buffer_rho(jbuf).ne.0.0e0) then
!              compute primitive contributions
               do mu=1, ngto
                  d(mu)=fx(i,mu)*fy(j,mu)*fz(k,mu)
                  dx(mu)=dfx(i,mu)*fy(j,mu)*fz(k,mu)
                  dy(mu)=fx(i,mu)*dfy(j,mu)*fz(k,mu)
                  dz(mu)=fx(i,mu)*fy(j,mu)*dfz(k,mu)
                  dx2(mu)=dfx2(i,mu)*fy(j,mu)*fz(k,mu)
                  dy2(mu)=fx(i,mu)*dfy2(j,mu)*fz(k,mu)
                  dz2(mu)=fx(i,mu)*fy(j,mu)*dfz2(k,mu)
               enddo
!              ****
!              transform to molecular orbitals
!              ****
               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
                  do mu=1,ngto
                     phi(imol)=phi(imol)+c(imol,mu)*d(mu)
                     dphi(imol,1)=dphi(imol,1)+c(imol,mu)*dx(mu)
                     dphi(imol,2)=dphi(imol,2)+c(imol,mu)*dy(mu)
                     dphi(imol,3)=dphi(imol,3)+c(imol,mu)*dz(mu)
                     d2phi(imol,1)=d2phi(imol,1)+c(imol,mu)*dx2(mu)
                     d2phi(imol,2)=d2phi(imol,2)+c(imol,mu)*dy2(mu)
                     d2phi(imol,3)=d2phi(imol,3)+c(imol,mu)*dz2(mu)
                  enddo
               enddo
               lmol=0
               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)
                     lapl=2.0d0*(dphi(imol,1)*dphi(jmol,1)+dphi(imol,2)*
     .                     dphi(jmol,2)+dphi(imol,3)*dphi(jmol,3))
                     laporb(lmol)=phi(imol)*(d2phi(jmol,1)+d2phi(jmol,2)
     .                          +d2phi(jmol,3))+phi(jmol)*(d2phi(imol,1)
     .                          +d2phi(imol,2)+d2phi(imol,3))+lapl
                     if(imol.ne.jmol) then
                         rhorb(lmol)=rhorb(lmol)+rhorb(lmol)
                         laporb(lmol)=laporb(lmol)+laporb(lmol)
                     endif
                  enddo
               enddo
!***
!              loop over molecular orbitals (one-int)
!              compute density and density laplacian
               rha=1.0d-20
               lapla=0.0d0
               gradfa=0.0d0
               do l=1,na
                  lmol=ia(l)
                  rha=rha+pa(l)*rhorb(lmol)
                  lapla=lapla+pa(l)*laporb(lmol)
               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))
               enddo
               if(nb.eq.0) then
                  rhb=rha
                  laplb=lapla
                  gradfb=gradfa
               else
!*    compute beta contribution (open shell case)
               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)
                     lapl=2.0d0*(dphi(imol,1)*dphi(jmol,1)+dphi(imol,2)*
     .                           dphi(jmol,2)+dphi(imol,3)*dphi(jmol,3))
                     laporb(lmol)=phi(imol)*(d2phi(jmol,1)+d2phi(jmol,2)
     .                          +d2phi(jmol,3))+phi(jmol)*(d2phi(imol,1)
     .                           +d2phi(imol,2)+d2phi(imol,3))+lapl    
                     if(imol.ne.jmol) then
                         rhorb(lmol)=rhorb(lmol)+rhorb(lmol)
                         laporb(lmol)=laporb(lmol)+laporb(lmol)
                     endif
                  enddo
               enddo
               rhb=1.0d-20
               laplb=0.0d0
               do l=1,nb
                  lmol=ib(l)
                  rhb=rhb+pb(l)*rhorb(lmol)
                  laplb=laplb+pb(l)*laporb(lmol)
               enddo
               gradfb=0.0d0
               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))
               enddo
               endif
               rh=rha+rhb
               lapl=-lapla-laplb
               gradf=(gradfa+gradfb)+tol*rh
               ratio=gradf*gradf*fermi*rh**(-16.0d0/3.0d0)
               elfd=1.d0/(1.0d0+ratio)
               buffer_elf(ibuf)=sngl(elfd)
               buffer_rho(ibuf)=sngl(rh)
               buffer_lap(ibuf)=sngl(lapl)
               endif
            enddo
         enddo
      enddo
      return
      end
