diff --git a/GoXC b/GoXC
index e91cb5a..52f4e57 100755
--- a/GoXC
+++ b/GoXC
@@ -13,6 +13,6 @@ then
   cp examples/basis."$1"."$2" input/basis
   cp examples/basis."$1"."$2" input/weight
   ./bin/IntPak
-  ./bin/xcDFT
+  ./bin/eDFT
 fi
 
diff --git a/examples/basis.Be.STO-3G b/examples/basis.Be.STO-3G
new file mode 100644
index 0000000..de39473
--- /dev/null
+++ b/examples/basis.Be.STO-3G
@@ -0,0 +1,13 @@
+1     3
+S   3   1.00
+      0.3016787069D+02       0.1543289673D+00
+      0.5495115306D+01       0.5353281423D+00
+      0.1487192653D+01       0.4446345422D+00
+S    3   1.00
+      0.1314833110D+01      -0.9996722919D-01       
+      0.3055389383D+00       0.3995128261D+00       
+      0.9937074560D-01       0.7001154689D+00       
+P    3   1.00
+      0.1314833110D+01      0.1559162750D+00
+      0.3055389383D+00      0.6076837186D+00
+      0.9937074560D-01      0.3919573931D+00
diff --git a/examples/molecule.H2 b/examples/molecule.H2
index a38d8a2..39795fb 100644
--- a/examples/molecule.H2
+++ b/examples/molecule.H2
@@ -1,5 +1,5 @@
-# nAt nEl nEla nElb nCore nRyd
-  2   2   1     1   0     0
+# nAt nEla nElb nCore nRyd
+  2   1     1   0     0
 # Znuc   x            y           z
   H     0.           0.         -0.7
   H     0.           0.         +0.7
diff --git a/input/basis b/input/basis
index b9ca7b5..eb58543 100644
--- a/input/basis
+++ b/input/basis
@@ -1,9 +1,29 @@
-1   3
-S   3   1.00
-     38.3600000              0.0238090        
-      5.7700000              0.1548910        
-      1.2400000              0.4699870        
+1   6
+S   8   1.00
+  17880.0000000              0.0007380        
+   2683.0000000              0.0056770        
+    611.5000000              0.0288830        
+    173.5000000              0.1085400        
+     56.6400000              0.2909070        
+     20.4200000              0.4483240        
+      7.8100000              0.2580260        
+      1.6530000              0.0150630        
+S   8   1.00
+  17880.0000000             -0.0001720        
+   2683.0000000             -0.0013570        
+    611.5000000             -0.0067370        
+    173.5000000             -0.0276630        
+     56.6400000             -0.0762080        
+     20.4200000             -0.1752270        
+      7.8100000             -0.1070380        
+      1.6530000              0.5670500        
 S   1   1.00
-      0.2976000              1.0000000        
+      0.4869000              1.0000000        
+P   3   1.00
+     28.3900000              0.0460870        
+      6.2700000              0.2401810        
+      1.6950000              0.5087440        
 P   1   1.00
-      1.2750000              1.0000000 
+      0.4317000              1.0000000        
+D   1   1.00
+      2.2020000              1.0000000  
diff --git a/input/dft b/input/dft
index ca9be70..90787c0 100644
--- a/input/dft
+++ b/input/dft
@@ -1,9 +1,9 @@
 # exchange rung: Hartree = 0, LDA = 1 (S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666
-  1 S51
+  666 HF
 # correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666
-  1 VWN5
+  0 HF
 # quadrature grid SG-n
-  3
+  1
 # Number of states in ensemble (nEns)
   1
 # Ensemble weights: wEns(1),...,wEns(nEns-1)
diff --git a/input/methods b/input/methods
index 8204ed8..341e057 100644
--- a/input/methods
+++ b/input/methods
@@ -5,7 +5,7 @@
 # CCD CCSD CCSD(T)
   F   F    F
 # CIS TDHF ADC
-  F   T    F
+  F   F    F
 # GF2 GF3
   F   F
 # G0W0 evGW qsGW
diff --git a/input/molecule b/input/molecule
index c78e87e..edeba31 100644
--- a/input/molecule
+++ b/input/molecule
@@ -1,4 +1,4 @@
 # nAt nEla nElb nCore nRyd
-  1   1     1   0     0
+  1   5     5   0     0
 # Znuc   x            y           z
-  He      0.0          0.0         0.0
+  Ne    0.0          0.0         0.0
diff --git a/input/options b/input/options
index 86d7547..d3d1e53 100644
--- a/input/options
+++ b/input/options
@@ -9,6 +9,6 @@
 # GF: maxSCF thresh  DIIS n_diis renormalization
       64     0.00001 T    10      3 
 # GW: maxSCF thresh  DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize
-      64     0.00001  T    20    F      F     T   F   F   F   F
+      64     0.00001  T    20    F      F     F   F   F   F   F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T
diff --git a/input/weight b/input/weight
index b9ca7b5..eb58543 100644
--- a/input/weight
+++ b/input/weight
@@ -1,9 +1,29 @@
-1   3
-S   3   1.00
-     38.3600000              0.0238090        
-      5.7700000              0.1548910        
-      1.2400000              0.4699870        
+1   6
+S   8   1.00
+  17880.0000000              0.0007380        
+   2683.0000000              0.0056770        
+    611.5000000              0.0288830        
+    173.5000000              0.1085400        
+     56.6400000              0.2909070        
+     20.4200000              0.4483240        
+      7.8100000              0.2580260        
+      1.6530000              0.0150630        
+S   8   1.00
+  17880.0000000             -0.0001720        
+   2683.0000000             -0.0013570        
+    611.5000000             -0.0067370        
+    173.5000000             -0.0276630        
+     56.6400000             -0.0762080        
+     20.4200000             -0.1752270        
+      7.8100000             -0.1070380        
+      1.6530000              0.5670500        
 S   1   1.00
-      0.2976000              1.0000000        
+      0.4869000              1.0000000        
+P   3   1.00
+     28.3900000              0.0460870        
+      6.2700000              0.2401810        
+      1.6950000              0.5087440        
 P   1   1.00
-      1.2750000              1.0000000 
+      0.4317000              1.0000000        
+D   1   1.00
+      2.2020000              1.0000000  
diff --git a/src/QuAcK/AO_values_grid.f90 b/src/QuAcK/AO_values_grid.f90
new file mode 100644
index 0000000..a5a72a6
--- /dev/null
+++ b/src/QuAcK/AO_values_grid.f90
@@ -0,0 +1,101 @@
+subroutine AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & 
+                          nGrid,root,AO,dAO)
+
+! Compute values of the AOs and their derivatives with respect to the cartesian coordinates
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas,nShell
+  double precision,intent(in)   :: CenterShell(maxShell,3)
+  integer,intent(in)            :: TotAngMomShell(maxShell)
+  integer,intent(in)            :: KShell(maxShell)
+  double precision,intent(in)   :: DShell(maxShell,maxK)
+  double precision,intent(in)   :: ExpShell(maxShell,maxK)
+  double precision,intent(in)   :: root(3,nGrid)
+  integer,intent(in)            :: nGrid
+
+! Local variables
+
+  integer                       :: atot,nShellFunction,a(3)
+  integer,allocatable           :: ShellFunction(:,:)
+  double precision              :: rASq,xA,yA,zA,norm_coeff,prim
+
+  integer                       :: iSh,iShF,iK,iG,iBas
+
+! Output variables
+
+  double precision,intent(out)  :: AO(nBas,nGrid)
+  double precision,intent(out)  :: dAO(3,nBas,nGrid)
+
+! Initialization
+
+  iBas = 0
+  AO(:,:)    = 0d0
+  dAO(:,:,:) = 0d0
+
+!------------------------------------------------------------------------
+! Loops over shells
+!------------------------------------------------------------------------
+  do iSh=1,nShell
+
+    atot = TotAngMomShell(iSh)
+    nShellFunction = (atot*atot + 3*atot + 2)/2
+    allocate(ShellFunction(1:nShellFunction,1:3))
+    call generate_shell(atot,nShellFunction,ShellFunction)
+
+    do iShF=1,nShellFunction
+
+      iBas = iBas + 1
+      a(:) = ShellFunction(iShF,:)
+
+      do iG=1,nGrid
+
+        xA = root(1,iG) - CenterShell(iSh,1)
+        yA = root(2,iG) - CenterShell(iSh,2)
+        zA = root(3,iG) - CenterShell(iSh,3)
+
+!       Calculate distance for exponential
+
+        rASq = xA**2 + yA**2 + zA**2
+
+!------------------------------------------------------------------------
+!         Loops over contraction degrees
+!-------------------------------------------------------------------------
+        do iK=1,KShell(iSh)
+          
+!         Calculate the exponential part
+
+          prim = DShell(iSh,iK)*norm_coeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq)        
+          AO(iBas,iG) = AO(iBas,iG) + prim
+
+          prim = -2d0*ExpShell(iSh,iK)*prim
+          dAO(:,iBas,iG) = dAO(:,iBas,iG) + prim
+        
+        enddo
+
+        dAO(1,iBas,iG) = xA**(a(1)+1)*yA**a(2)*zA**a(3)*dAO(1,iBas,iG)
+        if(a(1) > 0) dAO(1,iBas,iG) = dAO(1,iBas,iG) + dble(a(1))*xA**(a(1)-1)*yA**a(2)*zA**a(3)*AO(iBas,iG)
+        
+        dAO(2,iBas,iG) = xA**a(1)*yA**(a(2)+1)*zA**a(3)*dAO(2,iBas,iG)
+        if(a(2) > 0) dAO(2,iBas,iG) = dAO(2,iBas,iG) + dble(a(2))*xA**a(1)*yA**(a(2)-1)*zA**a(3)*AO(iBas,iG)
+        
+        dAO(3,iBas,iG) = xA**a(1)*yA**a(2)*zA**(a(3)+1)*dAO(3,iBas,iG)
+        if(a(3) > 0) dAO(3,iBas,iG) = dAO(3,iBas,iG) + dble(a(3))*xA**a(1)*yA**a(2)*zA**(a(3)-1)*AO(iBas,iG)
+
+!       Calculate polynmial part
+
+        AO(iBas,iG) = xA**a(1)*yA**a(2)*zA**a(3)*AO(iBas,iG)
+        
+      enddo
+      
+    enddo
+    deallocate(ShellFunction)
+  enddo
+!------------------------------------------------------------------------
+! End loops over shells 
+!------------------------------------------------------------------------
+
+end subroutine AO_values_grid
diff --git a/src/QuAcK/MO_values_grid.f90 b/src/QuAcK/MO_values_grid.f90
new file mode 100644
index 0000000..a0eeb35
--- /dev/null
+++ b/src/QuAcK/MO_values_grid.f90
@@ -0,0 +1,44 @@
+subroutine MO_values_grid(nBas,nGrid,c,AO,dAO,MO,dMO)
+
+! Compute values of the MOs and their derivatives with respect to the cartesian coordinates
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: c(nBas,nBas)
+  double precision,intent(in)   :: AO(nBas,nGrid)
+  double precision,intent(in)   :: dAO(3,nBas,nGrid)
+
+! Local variables
+
+  integer                       :: p,mu,iG
+
+! Output variables
+
+  double precision,intent(out)  :: MO(nBas,nGrid)
+  double precision,intent(out)  :: dMO(3,nBas,nGrid)
+
+! Initialization
+
+  MO(:,:)    = 0d0
+  dMO(:,:,:) = 0d0
+
+  do p=1,nBas
+    do mu=1,nBas
+      do iG=1,ngrid
+
+        MO(p,iG) = MO(p,iG) + c(mu,p)*AO(mu,iG)
+
+        dMO(1,p,iG) = dMO(1,p,iG) + c(mu,p)*dAO(1,mu,iG)
+        dMO(2,p,iG) = dMO(2,p,iG) + c(mu,p)*dAO(2,mu,iG)
+        dMO(3,p,iG) = dMO(3,p,iG) + c(mu,p)*dAO(3,mu,iG)
+
+      end do
+    end do
+  end do
+
+end subroutine MO_values_grid
diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90
index 01e0bae..2db13b5 100644
--- a/src/QuAcK/QuAcK.f90
+++ b/src/QuAcK/QuAcK.f90
@@ -76,6 +76,20 @@ program QuAcK
   double precision              :: dt
   logical                       :: doDrift
 
+  integer                       :: SGn
+  integer                       :: nRad
+  integer                       :: nAng
+  integer                       :: nGrid
+  double precision,allocatable  :: root(:,:)
+  double precision,allocatable  :: weight(:)
+  double precision,allocatable  :: AO(:,:)
+  double precision,allocatable  :: dAO(:,:,:)
+  double precision,allocatable  :: MO(:,:)
+  double precision,allocatable  :: dMO(:,:,:)
+  double precision,allocatable  :: rho(:)
+  double precision,allocatable  :: f(:)
+  double precision,allocatable  :: mu(:)
+
 ! Hello World
 
   write(*,*)
@@ -583,6 +597,36 @@ program QuAcK
 
   endif
 
+!------------------------------------------------------------------------
+! Basis set correction
+!------------------------------------------------------------------------
+
+!------------------------------------------------------------------------
+! Construct quadrature grid
+!------------------------------------------------------------------------
+
+  SGn = 1
+
+  call read_grid(SGn,nRad,nAng,nGrid)
+
+  allocate(root(ncart,nGrid),weight(nGrid))
+
+  call quadrature_grid(nRad,nAng,nGrid,root,weight)
+ 
+!------------------------------------------------------------------------
+! Calculate AO values at grid points
+!------------------------------------------------------------------------
+
+  allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid),MO(nBas,nGrid),dMO(ncart,nBas,nGrid))
+  allocate(rho(nGrid),f(nGrid),mu(nGrid))
+
+  call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO)
+  call density(nGrid,nBas,PHF(:,:,1),AO(:,:),rho(:))
+  call MO_values_grid(nBas,nGrid,cHF(:,:,1),AO,dAO,MO,dMO)
+  call f_grid(nBas,nO(1),nGrid,MO,ERI_MO_basis,f)
+  call mu_grid(nGrid,rho,f,mu)
+  call ec_srlda(nGrid,weight,rho,mu)
+
 !------------------------------------------------------------------------
 ! End of QuAcK
 !------------------------------------------------------------------------
diff --git a/src/QuAcK/density.f90 b/src/QuAcK/density.f90
index 90a17c1..e6bfd22 100644
--- a/src/QuAcK/density.f90
+++ b/src/QuAcK/density.f90
@@ -1,51 +1,38 @@
-subroutine density(doDrift,nBas,nWalk,P,gAO,dgAO,g,dg)
+subroutine density(nGrid,nBas,P,AO,rho)
 
-! Calculate the Green functions
+! Calculate one-electron density
 
   implicit none
+  include 'parameters.h'
 
 ! Input variables
 
-  logical,intent(in)            :: doDrift
-  integer,intent(in)            :: nBas,nWalk
-  double precision,intent(in)   :: P(nBas,nBas),gAO(nWalk,2,nBas),dgAO(nWalk,2,3,nBas)
+  double precision,parameter    :: thresh = 1d-15
+
+  integer,intent(in)            :: nGrid
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: P(nBas,nBas)
+  double precision,intent(in)   :: AO(nBas,nGrid)
 
 ! Local variables
 
-  integer                       :: iW,iEl,ixyz,mu,nu
+  integer                       :: iG,mu,nu
 
 ! Output variables
 
-  double precision,intent(out)  :: g(nWalk,2),dg(nWalk,2,3)
+  double precision,intent(out)  :: rho(nGrid) 
 
-  g = 0d0
-  do iW=1,nWalk
-   do iEl=1,2
-     do mu=1,nBas
-       do nu=1,nBas
-         g(iW,iEl) = g(iW,iEl) + gAO(iW,iEl,mu)*P(mu,nu)*gAO(iW,iEl,nu)
-       enddo
-     enddo
-   enddo
-  enddo
-
-  if(doDrift) then
-
-    dg = 0d0
-    do iW=1,nWalk
-      do iEl=1,2
-        do ixyz=1,3
-          do mu=1,nBas
-            do nu=1,nBas
-              dg(iW,iEl,ixyz) = dg(iW,iEl,ixyz)                                & 
-                              + P(mu,nu)*(dgAO(iW,iEl,ixyz,mu)*gAO(iW,iEl,nu)  &
-                              +            gAO(iW,iEl,mu)*dgAO(iW,iEl,ixyz,nu))
-            enddo
-          enddo
-        enddo
+  rho(:) = 0d0
+  do iG=1,nGrid
+    do mu=1,nBas
+      do nu=1,nBas
+        rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG)
       enddo
     enddo
+  enddo
 
-  endif
+! do iG=1,nGrid
+!   rho(iG) = max(rho(iG),thresh)
+! enddo
 
 end subroutine density
diff --git a/src/eDFT/dft_grid.f90 b/src/QuAcK/dft_grid.f
similarity index 100%
rename from src/eDFT/dft_grid.f90
rename to src/QuAcK/dft_grid.f
diff --git a/src/QuAcK/ec_srlda.f90 b/src/QuAcK/ec_srlda.f90
new file mode 100644
index 0000000..85b653a
--- /dev/null
+++ b/src/QuAcK/ec_srlda.f90
@@ -0,0 +1,46 @@
+subroutine ec_srlda(nGrid,weight,rho,mu)
+
+! Compute sr-lda ec
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rho(nGrid)
+  double precision,intent(in)   :: mu(nGrid)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: r
+  double precision              :: rs
+  double precision              :: ecsr
+  double precision              :: ec
+  double precision,parameter    :: thres = 1d-15
+
+! Initialization
+
+  ec = 0d0
+
+  do iG=1,ngrid
+
+    r = max(0d0,rho(iG))
+
+    if(r > thres) then
+
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+
+      call srlda(rs,mu(iG),ecsr)
+
+      ec = ec + weight(iG)*ecsr*rho(iG)
+
+    end if
+
+  end do
+
+  print*, 'ec = ',ec
+
+end subroutine ec_srlda
diff --git a/src/QuAcK/excitation_density_SOSEX_RI.f90 b/src/QuAcK/excitation_density_SOSEX_RI.f90
new file mode 100644
index 0000000..8f3f6b5
--- /dev/null
+++ b/src/QuAcK/excitation_density_SOSEX_RI.f90
@@ -0,0 +1,65 @@
+subroutine excitation_density_SOSEX_RI(nBas,nC,nO,nR,nS,c,G,XpY,rho)
+
+! Compute excitation densities
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas,nC,nO,nR,nS
+  double precision,intent(in)   :: c(nBas,nBas),G(nBas,nBas,nBas,nBas),XpY(nS,nS)
+
+! Local variables
+
+  double precision,allocatable  :: scr(:,:,:)
+  integer                       :: mu,nu,la,si,ia,jb,x,y,j,b
+
+! Output variables
+
+  double precision,intent(out)  :: rho(nBas,nBas,nS)
+
+! Memory allocation
+  allocate(scr(nBas,nBas,nS))
+
+  rho(:,:,:) = 0d0
+  do nu=1,nBas
+    do si=1,nBas
+      do ia=1,nS
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+            rho(nu,si,ia) = rho(nu,si,ia) + c(nu,j)*XpY(ia,jb)*c(si,b)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+ 
+  scr(:,:,:) = 0d0   
+  do mu=1,nBas
+    do la=1,nBas
+      do ia=1,nS
+        do nu=1,nBas
+          do si=1,nBas
+            scr(mu,la,ia) = scr(mu,la,ia) + G(mu,nu,la,si)*rho(nu,si,ia)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  rho(:,:,:) = 0d0   
+  do ia=1,nS
+    do x=nC+1,nBas-nR
+      do y=nC+1,nBas-nR
+        do mu=1,nBas
+          do la=1,nBas
+            rho(x,y,ia) = rho(x,y,ia) + c(mu,x)*scr(mu,la,ia)*c(la,y)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+end subroutine excitation_density_SOSEX_RI
diff --git a/src/QuAcK/f_grid.f90 b/src/QuAcK/f_grid.f90
new file mode 100644
index 0000000..63bd3f8
--- /dev/null
+++ b/src/QuAcK/f_grid.f90
@@ -0,0 +1,44 @@
+subroutine f_grid(nBas,nO,nGrid,MO,ERI,f)
+
+! Compute f
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: MO(nBas,nGrid)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer                       :: p,q
+  integer                       :: i,j
+  integer                       :: iG
+
+! Output variables
+
+  double precision,intent(out)  :: f(nGrid)
+
+! Initialization
+
+  f(:) = 0d0
+
+  do p=1,nBas
+    do q=1,nBas
+      do i=1,nO
+        do j=1,nO
+          do iG=1,ngrid
+
+            f(iG) = f(iG) + MO(p,iG)*MO(q,iG)*ERI(p,q,i,j)*MO(i,iG)*MO(j,iG)
+
+          end do
+        end do
+      end do
+    end do
+  end do
+
+end subroutine f_grid
diff --git a/src/QuAcK/mu_grid.f90 b/src/QuAcK/mu_grid.f90
new file mode 100644
index 0000000..23a6ff2
--- /dev/null
+++ b/src/QuAcK/mu_grid.f90
@@ -0,0 +1,46 @@
+subroutine mu_grid(nGrid,rho,f,mu)
+
+! Compute mu
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: rho(nGrid)
+  double precision,intent(in)   :: f(nGrid)
+
+! Local variables
+
+  integer                       :: iG
+  double precision,parameter    :: thres = 1d-10
+  double precision              :: n2
+
+! Output variables
+
+  double precision,intent(out)  :: mu(nGrid)
+
+! Initialization
+
+  mu(:) = 0d0
+
+  do iG=1,ngrid
+
+    n2 = rho(iG)**2
+
+    if(abs(n2) > thres) then
+
+      mu(iG) =  f(iG)/n2
+   
+    else
+    
+      mu(iG) = 1d0/thres
+
+    end if
+
+  end do
+
+  mu(:) = 0.5d0*sqrt(pi)*mu(:)
+
+end subroutine mu_grid
diff --git a/src/QuAcK/quadrature_grid.f90 b/src/QuAcK/quadrature_grid.f90
new file mode 100644
index 0000000..420e80a
--- /dev/null
+++ b/src/QuAcK/quadrature_grid.f90
@@ -0,0 +1,77 @@
+subroutine quadrature_grid(nRad,nAng,nGrid,root,weight)
+
+! Build roots and weights of quadrature grid
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)                 :: nRad,nAng,nGrid
+
+! Local variables
+
+  integer                            :: i,j,k
+  double precision                   :: scale
+  double precision,allocatable       :: Radius(:)
+  double precision,allocatable       :: RadWeight(:)
+  double precision,allocatable       :: XYZ(:,:)
+  double precision,allocatable       :: XYZWeight(:)
+
+! Output variables
+
+  double precision,intent(out)       :: root(3,nGrid)
+  double precision,intent(out)       :: weight(nGrid)
+
+! Memory allocation
+
+  allocate(Radius(nRad),RadWeight(nRad),XYZ(3,nAng),XYZWeight(nAng))
+
+! Findthe radial grid
+
+  scale = 1d0
+  call EulMac(Radius,RadWeight,nRad,scale)
+
+  write(*,20)
+  write(*,30)
+  write(*,20)
+  do i = 1,nRad
+     write(*,40) i,Radius(i),RadWeight(i)
+  end do
+  write(*,20)
+  write(*,*)
+
+! Find the angular grid
+
+  call Lebdev(XYZ,XYZWeight,nAng)
+
+  write(*,20)
+  write(*,50)
+  write(*,20)
+  do j = 1,nAng
+     write(*,60) j,(XYZ(k,j),k=1,3), XYZWeight(j)
+  end do
+  write(*,20)
+
+! Form the roots and weights
+
+  k = 0
+  do i=1,nRad
+    do j=1,nAng
+       k = k + 1
+       root(:,k) = Radius(i)*XYZ(:,j)
+       weight(k) = RadWeight(i)*XYZWeight(j)
+    enddo
+  enddo
+ 
+! Compute values of the basis functions (and the its gradient if required) at each grid point
+
+20 format(T2,58('-'))
+30 format(T20,'Radial Quadrature',/, &
+          T6,'I',T26,'Radius',T50,'Weight')
+40 format(T3,I4,T18,F17.10,T35,F25.10)
+50 format(T20,'Angular Quadrature',/, &
+          T6,'I',T19,'X',T29,'Y',T39,'Z',T54,'Weight')
+60 format(T3,I4,T13,3F10.5,T50,F10.5)
+
+end subroutine quadrature_grid
diff --git a/src/QuAcK/read_grid.f90 b/src/QuAcK/read_grid.f90
new file mode 100644
index 0000000..5cf4391
--- /dev/null
+++ b/src/QuAcK/read_grid.f90
@@ -0,0 +1,47 @@
+subroutine read_grid(SGn,nRad,nAng,nGrid)
+
+! Read grid type
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: SGn
+
+! Output variables
+
+  integer,intent(out)           :: nRad
+  integer,intent(out)           :: nAng
+  integer,intent(out)           :: nGrid
+
+  write(*,*)'----------------------------------------------------------'
+  write(*,'(A22,I1)')'  Quadrature grid: SG-',SGn
+  write(*,*)'----------------------------------------------------------'
+
+  select case (SGn)
+
+    case(0)
+     nRad = 23
+     nAng = 170
+
+    case(1)
+      nRad = 50
+      nAng = 194
+
+    case(2)
+      nRad = 75
+      nAng = 302
+
+    case(3)
+      nRad = 99
+      nAng = 590
+
+    case default
+      call print_warning('!!! Quadrature grid not available !!!')
+      stop
+
+  end select
+
+  nGrid = nRad*nAng
+
+end subroutine read_grid
diff --git a/src/QuAcK/srlda.f90 b/src/QuAcK/srlda.f90
new file mode 100644
index 0000000..aa8609e
--- /dev/null
+++ b/src/QuAcK/srlda.f90
@@ -0,0 +1,125 @@
+subroutine srlda(rs,mu,ecsr)
+
+! Correlation energy of a spin unpolarized uniform electron gas 
+! with short-range interaction erfc(mu*r)/r
+! See Zecca et al. PRB 70, 205127 (2004)
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)  :: rs
+  double precision,intent(in)  :: mu
+
+! Local variables
+
+  double precision             :: ec
+  double precision             :: cf
+  double precision             :: b1
+  double precision             :: b2
+  double precision             :: b3
+  double precision             :: b4
+  double precision             :: a0
+  double precision             :: bb
+  double precision             :: m1
+
+! Ouput variables
+
+  double precision,intent(out) :: ecsr
+
+! Compute PW LDA correlation energy
+
+  call ecPW(rs,0d0,ec)
+
+! Define various stuff
+
+  cf = (9d0*pi/4d0)**(1d0/3d0)
+  bb = 1.27329d0    
+  m1 = 0.0357866d0  
+  a0 = ec
+  b3 = bb*rs**(7d0/2d0)
+  b2 = -3d0/2d0/pi/cf*rs/a0
+  b1 = (b3-1d0/sqrt(3d0*pi)*rs**(3d0/2d0)/a0)/b2
+  b4 = -a0*b1*rs**3/m1
+
+! Compute short-range correlation energy
+
+  ecsr = a0*(1d0 + b1*mu)/(1d0 + b1*mu+b2*mu**2 + b3*mu**3 + b4*mu**4)
+
+end subroutine srlda
+
+!==========================================================================================
+
+subroutine ecPW(x,y,ec)
+
+! Correlation energy of the 3D electron gas of density rs and spin polarization z
+! Perdew & Wang, PRB 45, 13244 (1992)
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables 
+
+  double precision,intent(in)  :: x
+  double precision,intent(in)  :: y
+
+! Local variables 
+
+  double precision             :: f02
+  double precision             :: ff
+  double precision             :: aaa
+  double precision             :: G
+  double precision             :: ec0
+  double precision             :: ec1
+  double precision             :: alfac
+  
+! Ouput variables
+
+  double precision,intent(out) :: ec
+
+  f02 = 4d0/(9d0*(2d0**(1d0/3d0) - 1d0))
+
+  ff = ((1d0+y)**(4d0/3d0) + (1d0-y)**(4d0/3d0)-2d0)/(2d0**(4d0/3d0) - 2d0)
+
+  aaa = (1d0 - log(2d0))/pi**2
+
+  call GPW(x,aaa,0.21370d0,7.5957d0,3.5876d0,1.6382d0,0.49294d0,G)
+  ec0 = G
+
+  aaa=aaa/2d0
+  call GPW(x,aaa,0.20548d0,14.1189d0,6.1977d0,3.3662d0,0.62517d0,G)
+  ec1 = G
+
+  call GPW(x,0.016887d0,0.11125d0,10.357d0,3.6231d0,0.88026d0,0.49671d0,G)
+  alfac = -G
+
+  ec = ec0 + alfac*ff/f02*(1d0 - y**4) + (ec1 - ec0)*ff*y**4
+
+end subroutine ecPW
+
+!==========================================================================================
+
+subroutine GPW(x,Ac,alfa1,beta1,beta2,beta3,beta4,G)
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)  :: Ac
+  double precision,intent(in)  :: alfa1
+  double precision,intent(in)  :: beta1
+  double precision,intent(in)  :: beta2
+  double precision,intent(in)  :: beta3
+  double precision,intent(in)  :: beta4
+  double precision,intent(in)  :: x
+
+! Ouput variables
+
+  double precision,intent(out) :: G
+
+  G = -2d0*Ac*(1d0 + alfa1*x)*log(1d0 & 
+    + 1d0/(2d0*Ac*(beta1*sqrt(x) + beta2*x + beta3*x*sqrt(x) + beta4*x**2)))
+
+end subroutine GPW
diff --git a/src/eDFT/DIIS_extrapolation.f90 b/src/eDFT/DIIS_extrapolation.f90
deleted file mode 100644
index e53db06..0000000
--- a/src/eDFT/DIIS_extrapolation.f90
+++ /dev/null
@@ -1,54 +0,0 @@
-subroutine DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_inout)
-
-! Perform DIIS extrapolation
-
-  implicit none
-
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: n_err,n_e
-  double precision,intent(in)   :: error_in(n_err),error(n_err,n_diis),e(n_e,n_diis)
-
-! Local variables
-
-  double precision,allocatable  :: A(:,:),b(:),w(:)
-
-! Output variables
-
-  double precision,intent(out)  :: rcond
-  integer,intent(inout)         :: n_diis
-  double precision,intent(inout):: e_inout(n_e)
-
-! Memory allocaiton
-
-  allocate(A(n_diis+1,n_diis+1),b(n_diis+1),w(n_diis+1))
-
-! Update DIIS "history"
-
-  call prepend(n_err,n_diis,error,error_in)
-  call prepend(n_e,n_diis,e,e_inout)
-
-!  Build A matrix
-
-  A(1:n_diis,1:n_diis) = matmul(transpose(error),error)
-
-  A(1:n_diis,n_diis+1) = -1d0
-  A(n_diis+1,1:n_diis) = -1d0
-  A(n_diis+1,n_diis+1) = +0d0
-
-! Build x matrix
-
-  b(1:n_diis) = +0d0
-  b(n_diis+1) = -1d0
-
-! Solve linear system
-
-  call linear_solve(n_diis+1,A,b,w,rcond)
-
-! Extrapolate
-
-    e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis)))
-
-end subroutine DIIS_extrapolation
diff --git a/src/eDFT/Makefile b/src/eDFT/Makefile
index a72dc51..2c3316c 100644
--- a/src/eDFT/Makefile
+++ b/src/eDFT/Makefile
@@ -1,12 +1,17 @@
 IDIR =../../include
 BDIR =../../bin
 ODIR = obj
+OODIR = ../IntPak/obj
 SDIR =.
-FC = gfortran -I$(IDIR) 
+FC = gfortran -I$(IDIR) -Wall -g -Wno-unused -Wno-unused-dummy-argument
 ifeq ($(DEBUG),1)
 FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
 else
-FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O2 
+FFLAGS = -O3
+endif
+
+ifeq ($(PROFIL),1)
+  FC += -p -fno-inline
 endif
 
 LIBS = ~/Dropbox/quack/lib/*.a
@@ -18,7 +23,7 @@ SRC = $(wildcard *.f)
 
 OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC))
 
-$(ODIR)/%.o: %.f90
+$(ODIR)/%.o: %.f90 
 	$(FC) -c -o $@ $< $(FFLAGS) 
 
 $(ODIR)/%.o: %.f
@@ -29,7 +34,9 @@ $(BDIR)/eDFT: $(OBJ)
 
 debug: 
 	DEBUG=1 make $(BDIR)/eDFT
-#DEBUG=1 make clean $(BDIR)/eDFT
+
+profil: 
+	PROFIL=1 make $(BDIR)/eDFT
 
 clean:
 	rm -f $(ODIR)/*.o $(BDIR)/eDFT $(BDIR)/debug
diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90
index 9db226b..9ebae91 100644
--- a/src/eDFT/eDFT.f90
+++ b/src/eDFT/eDFT.f90
@@ -2,9 +2,11 @@ program eDFT
 
 ! exchange-correlation density-functional theory calculations
 
+  implicit none
   include 'parameters.h'
 
-  integer                       :: nNuc,nBas,nEl(nspin),nO(nspin),nV(nspin)
+  integer                       :: nNuc,nBas
+  integer                       :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin)
   double precision              :: ENuc,EKS
 
   double precision,allocatable  :: ZNuc(:),rNuc(:,:)
@@ -35,7 +37,9 @@ program eDFT
 
   integer                       :: maxSCF,max_diis
   double precision              :: thresh
-  logical                       :: DIIS,guess_type,ortho_type
+  logical                       :: DIIS
+  integer                       :: guess_type
+  integer                       :: ortho_type
 
 ! Hello World
 
@@ -50,12 +54,14 @@ program eDFT
 !------------------------------------------------------------------------
 
 ! Read number of atoms, number of electrons of the system
+! nC   = number of core orbitals
 ! nO   = number of occupied orbitals
 ! nV   = number of virtual orbitals (see below)
+! nR   = number of Rydberg orbitals 
 ! nBas = number of basis functions (see below)
 !      = nO + nV
 
-  call read_molecule(nNuc,nEl,nO)
+  call read_molecule(nNuc,nEl(:),nO(:),nC(:),nR(:))
   allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
 
 ! Read geometry
@@ -71,23 +77,6 @@ program eDFT
 
   call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
 
-!------------------------------------------------------------------------
-! Read one- and two-electron integrals
-!------------------------------------------------------------------------
-
-! Memory allocation for one- and two-electron integrals
-
-  allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), &
-           ERI(nBas,nBas,nBas,nBas))
-
-!   Read integrals
-
-  call read_integrals(nBas,S,T,V,Hc,ERI)  
-
-! Orthogonalization X = S^(-1/2)
-
-  call orthogonalization_matrix(nBas,S,X)
-
 !------------------------------------------------------------------------
 ! DFT options
 !------------------------------------------------------------------------
@@ -95,9 +84,24 @@ program eDFT
 ! Allocate ensemble weights
 
   allocate(wEns(maxEns))
-
   call read_options(x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type)
 
+!------------------------------------------------------------------------
+! Read one- and two-electron integrals
+!------------------------------------------------------------------------
+
+! Memory allocation for one- and two-electron integrals
+
+  allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI(nBas,nBas,nBas,nBas))
+
+! Read integrals
+
+  call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI)  
+
+! Orthogonalization X = S^(-1/2)
+
+  call orthogonalization_matrix(ortho_type,nBas,S,X)
+
 !------------------------------------------------------------------------
 ! Construct quadrature grid
 !------------------------------------------------------------------------
@@ -111,8 +115,7 @@ program eDFT
 !------------------------------------------------------------------------
 
   allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid))
-  call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
-                      nGrid,root,AO,dAO)
+  call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO)
 
 !------------------------------------------------------------------------
 ! Compute KS energy
diff --git a/src/eDFT/orthogonalization_matrix.f90 b/src/eDFT/orthogonalization_matrix.f90
index ed837d4..15ea4ac 100644
--- a/src/eDFT/orthogonalization_matrix.f90
+++ b/src/eDFT/orthogonalization_matrix.f90
@@ -1,12 +1,12 @@
-subroutine orthogonalization_matrix(nBas,S,X)
+subroutine orthogonalization_matrix(ortho_type,nBas,S,X)
 
-! Compute the orthogonalization matrix X = S^(-1/2)
+! Compute the orthogonalization matrix X
 
   implicit none
 
 ! Input variables
 
-  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nBas,ortho_type
   double precision,intent(in)   :: S(nBas,nBas)
 
 ! Local variables
@@ -23,30 +23,87 @@ subroutine orthogonalization_matrix(nBas,S,X)
 
   debug = .false.
 
+! Type of orthogonalization ortho_type
+!
+!  1 = Lowdin
+!  2 = Canonical
+!  3 = SVD
+!
+
   allocate(Uvec(nBas,nBas),Uval(nBas))
 
-  write(*,*)
-  write(*,*) ' *** Lowdin orthogonalization X = S^(-1/2)  *** '
-  write(*,*)
+  if(ortho_type == 1) then
 
-  Uvec = S
-  call diagonalize_matrix(nBas,Uvec,Uval)
+    write(*,*)
+    write(*,*) ' Lowdin orthogonalization'
+    write(*,*)
 
-  do i=1,nBas
+    Uvec = S
+    call diagonalize_matrix(nBas,Uvec,Uval)
 
-    if(Uval(i) > thresh) then 
+    do i=1,nBas
 
-      Uval(i) = 1d0/sqrt(Uval(i))
+      if(Uval(i) > thresh) then 
 
-    else
+        Uval(i) = 1d0/sqrt(Uval(i))
 
-      write(*,*) 'Eigenvalue',i,'too small for Lowdin orthogonalization'
+      else
 
-    endif
+        write(*,*) 'Eigenvalue',i,'too small for Lowdin orthogonalization'
 
-  enddo
-  
-  call ADAt(nBas,Uvec,Uval,X)
+      endif
+
+    enddo
+    
+    call ADAt(nBas,Uvec,Uval,X)
+
+  elseif(ortho_type == 2) then
+
+    write(*,*)
+    write(*,*) 'Canonical orthogonalization'
+    write(*,*)
+
+    Uvec = S
+    call diagonalize_matrix(nBas,Uvec,Uval)
+
+    do i=1,nBas
+
+      if(Uval(i) > thresh) then 
+
+        Uval(i) = 1d0/sqrt(Uval(i))
+
+      else
+
+        write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization'
+
+      endif
+
+    enddo
+    
+    call AD(nBas,Uvec,Uval)
+    X = Uvec
+ 
+  elseif(ortho_type == 3) then
+
+    write(*,*)
+    write(*,*) ' SVD-based orthogonalization NYI'
+    write(*,*)
+   
+!   Uvec = S
+!   call diagonalize_matrix(nBas,Uvec,Uval)
+
+!   do i=1,nBas
+!     if(Uval(i) > thresh) then 
+!       Uval(i) = 1d0/sqrt(Uval(i))
+!     else
+!       write(*,*) 'Eigenvalue',i,'too small for canonical orthogonalization'
+!     endif
+!   enddo
+!   
+!   call AD(nBas,Uvec,Uval)
+!   X = Uvec
+ 
+  endif
 
 ! Print results
 
diff --git a/src/utils/Makefile b/src/utils/Makefile
index 09a8268..d30e9a9 100644
--- a/src/utils/Makefile
+++ b/src/utils/Makefile
@@ -5,7 +5,7 @@ ODIR = obj
 OODIR = ../IntPak/obj
 SDIR =.
 FC = gfortran -I$(IDIR) 
-AR = ar
+AR = libtool
 ifeq ($(DEBUG),1)
 FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
 else
@@ -28,7 +28,7 @@ $(ODIR)/%.o: %.f
 	$(FC) -c -o $@ $< $(FFLAGS) 
 
 $(LDIR)/utils.a: $(OBJ)
-	$(AR) -cru $@ $^ 
+	$(AR) -static -o $@ $^ 
 
 debug: 
 	DEBUG=1 make $(LDIR)/utils.a
diff --git a/src/utils/map_module.f90 b/src/utils/map_module.f90
deleted file mode 100644
index 98e7347..0000000
--- a/src/utils/map_module.f90
+++ /dev/null
@@ -1,902 +0,0 @@
-module map_module
-
-! A map is an array of maps (cache_maps)
-! A cache map is an array of keys and values sorted by keys
-! A cache map has its own OpenMP lock
-! To access a (key,value) pair in the map, the
-! index of the cache_map in the map array is obtained
-! by removing the first 15 bits of the key.
-! The key in the cache_map is composed of the first
-! 15 bits of the key. Therefore, it can be stored
-! as integer*2 and is found by applying the map_mask
-! to the initial key. The element are found in the
-! cache_map using a binary search
-!
-! When using the map_update subroutine to build the map,
-! the map_merge subroutine
-! should be called before getting data from the map.
-
- use omp_lib
-
- integer, parameter             :: integral_kind = 8
-
- integer, parameter             :: cache_key_kind = 2
- integer, parameter             :: cache_map_size_kind = 4
-
- integer, parameter             :: key_kind = 8
- integer, parameter             :: map_size_kind = 8
-
- integer,   parameter           :: map_shift = 15
- integer*8, parameter           :: map_mask  = ibset(0_8,15)-1_8
-
- type cache_map_type
-  real(integral_kind), pointer   :: value(:)
-  integer(cache_key_kind), pointer :: key(:)
-  logical                        :: sorted
-  integer(cache_map_size_kind)   :: map_size
-  integer(cache_map_size_kind)   :: n_elements
-  integer(omp_lock_kind)         :: lock
- end type cache_map_type
-
- type map_type
-  type(cache_map_type), allocatable :: map(:)
-  real(integral_kind), pointer   :: consolidated_value(:)
-  integer(cache_key_kind), pointer :: consolidated_key(:)
-  integer*8, pointer             :: consolidated_idx(:)
-  logical                        :: sorted
-  logical                        :: consolidated
-  integer(map_size_kind)         :: map_size
-  integer(map_size_kind)         :: n_elements
-  integer(omp_lock_kind)         :: lock
- end type map_type
-
-end module map_module
-
-
-double precision function map_mb(map)
-  use map_module
-  use omp_lib
-  implicit none
-  type (map_type), intent(in)    :: map
-  integer(map_size_kind)         :: i
-
-  map_mb = dble(8+map_size_kind+map_size_kind+omp_lock_kind+4)
-  do i=0,map%map_size
-    map_mb = map_mb + dble(map%map(i)%map_size*(cache_key_kind+integral_kind) +&
-        8+8+4+cache_map_size_kind+cache_map_size_kind+omp_lock_kind)
-  enddo
-  map_mb = map_mb / (1024.d0*1024.d0)
-end
-
-subroutine cache_map_init(map,sze)
-  use map_module
-  implicit none
-  type (cache_map_type), intent(inout) :: map
-  integer(cache_map_size_kind)   :: sze
-  call omp_set_lock(map%lock)
-  map%n_elements = 0_8
-  map%map_size = 0_8
-  map%sorted = .True.
-  NULLIFY(map%value, map%key)
-  call cache_map_reallocate(map,sze)
-  call omp_unset_lock(map%lock)
-end
-
-subroutine map_init(map,keymax)
-  use map_module
-  implicit none
-  integer*8, intent(in)          :: keymax
-  type (map_type), intent(inout) :: map
-  integer(map_size_kind)         :: i
-  integer(cache_map_size_kind)   :: sze
-  integer                        :: err
-
-  call omp_init_lock(map%lock)
-  call omp_set_lock(map%lock)
-
-  map%n_elements = 0_8
-  map%map_size = shiftr(keymax,map_shift)
-  map%consolidated = .False.
-
-  allocate(map%map(0_8:map%map_size),stat=err)
-  if (err /= 0) then
-    print *,  'Unable to allocate map'
-    stop 5
-  endif
-  sze = 2
-  do i=0_8,map%map_size
-    call omp_init_lock(map%map(i)%lock)
-  enddo
-  !$OMP PARALLEL DEFAULT(NONE) SHARED(map,sze) PRIVATE(i)
-  !$OMP DO SCHEDULE(STATIC,512)
-  do i=0_8,map%map_size
-    call cache_map_init(map%map(i),sze)
-  enddo
-  !$OMP ENDDO
-  !$OMP END PARALLEL
-  map%sorted = .True.
-
-  call omp_unset_lock(map%lock)
-
-end
-
-subroutine cache_map_reallocate(map,sze)
-  use map_module
-  implicit none
-  integer(cache_map_size_kind), intent(in) :: sze
-  type (cache_map_type), intent(inout) :: map
-
-  integer(cache_key_kind), pointer :: key_new(:)
-  real(integral_kind), pointer   :: value_new(:)
-  integer(map_size_kind)         :: i
-  integer                        :: err
-  !DIR$ ATTRIBUTES ALIGN : 64    :: key_new, value_new
-
-  if (sze < map%n_elements) then
-    print *,  'Unable to resize map : map too large'
-    stop 3
-  endif
-
-  ! Resize keys
-  allocate( key_new(sze), stat=err )
-  if (err /= 0) then
-    print *,  'Unable to allocate map', sze
-    stop 1
-  endif
-  if (associated(map%key)) then
-    do i=1_8,min(size(map%key),map%n_elements)
-      key_new(i) = map%key(i)
-    enddo
-    deallocate(map%key)
-  endif
-
-  ! Resize values
-  allocate( value_new(sze), stat=err )
-  if (err /= 0) then
-    print *,  'Unable to allocate map', sze
-    stop 2
-  endif
-  if (associated(map%value)) then
-    do i=1_8,min(size(map%key),map%n_elements)
-      value_new(i) = map%value(i)
-    enddo
-    deallocate(map%value)
-  endif
-
-  ! Set new pointers
-  map%key => key_new
-  map%value => value_new
-  map%map_size = sze
-
-end
-
-
-subroutine cache_map_deinit(map)
-  use map_module
-  implicit none
-  type (cache_map_type), intent(inout) :: map
-
-  integer                        :: err
-
-  if (associated( map % value )) then
-    deallocate( map % value, stat=err )
-    if (err /= 0) then
-      print *,  'Unable to deallocate map'
-      stop 2
-    endif
-    NULLIFY(map%value)
-  endif
-
-  if (associated( map % key )) then
-    deallocate( map % key, stat=err )
-    if (err /= 0) then
-      print *,  'Unable to deallocate map'
-      stop 4
-    endif
-    NULLIFY(map%key)
-  endif
-
-  map%n_elements = 0_8
-  map%map_size   = 0_8
-  call omp_destroy_lock(map%lock)
-end
-
-subroutine map_deinit(map)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  integer                        :: err
-  integer(map_size_kind)         :: i
-
-  if (allocated( map % map )) then
-    do i=0_8, map%map_size
-      call cache_map_deinit(map%map(i))
-    enddo
-    deallocate( map % map, stat=err )
-    if (err /= 0) then
-      print *,  'Unable to deallocate map'
-      stop 6
-    endif
-  endif
-
-  map%n_elements = 0_8
-  map%map_size   = 0_8
-  call omp_destroy_lock(map%lock)
-end
-
-subroutine cache_map_sort(map)
-  use map_module
-  implicit none
-  type (cache_map_type), intent(inout) :: map
-  integer(cache_map_size_kind), allocatable :: iorder(:)
-  integer(cache_map_size_kind)   :: i
-  !DIR$ ATTRIBUTES ALIGN : 64    :: iorder
-
-  if (.not.map%sorted) then
-    allocate(iorder(map%n_elements))
-    do i=1,map%n_elements
-      iorder(i) = i
-    enddo
-    if (cache_key_kind == 2) then
-      call i2radix_sort(map%key,iorder,map%n_elements,-1)
-    else if (cache_key_kind == 4) then
-      call iradix_sort(map%key,iorder,map%n_elements,-1)
-    else if (cache_key_kind == 8) then
-      call i8radix_sort(map%key,iorder,map%n_elements,-1)
-    endif
-    if (integral_kind == 4) then
-      call set_order(map%value,iorder,map%n_elements)
-    else if (integral_kind == 8) then
-      call dset_order(map%value,iorder,map%n_elements)
-    endif
-    deallocate(iorder)
-    map%sorted = .True.
-  endif
-
-end
-
-subroutine map_sort(map)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  integer(map_size_kind)         :: i
-
-  if (.not.map%sorted) then
-    !$OMP PARALLEL DO SCHEDULE(static,1024) DEFAULT(SHARED) PRIVATE(i)
-    do i=0_8,map%map_size
-      call omp_set_lock(map%map(i)%lock)
-      call cache_map_sort(map%map(i))
-      call omp_unset_lock(map%map(i)%lock)
-    enddo
-    !$OMP END PARALLEL  DO
-    map%sorted = .True.
-  endif
-
-end
-
-subroutine cache_map_merge(map)
-  use map_module
-  implicit none
-  type (cache_map_type), intent(inout) :: map
-  integer(cache_key_kind)        :: prev_key
-  integer(cache_map_size_kind)   :: i, j
-
-  call cache_map_sort(map)
-  prev_key = -1_8
-  j=0
-  do i=1,map%n_elements
-    if (map%key(i) /= prev_key) then
-      j = j+1
-      map%value(j) = map%value(i)
-      map%key(j) = map%key(i)
-      prev_key = map%key(i)
-    else
-      map%value(j) = map%value(j)+map%value(i)
-    endif
-  enddo
-  map%n_elements = j
-
-end
-
-subroutine cache_map_unique(map)
-  use map_module
-  implicit none
-  type (cache_map_type), intent(inout) :: map
-  integer(cache_key_kind)        :: prev_key
-  integer(cache_map_size_kind)   :: i, j
-
-  call cache_map_sort(map)
-  prev_key = -1_8
-  j=0
-  do i=1,map%n_elements
-    if (map%key(i) /= prev_key) then
-      j = j+1
-      map%value(j) = map%value(i)
-      map%key(j) = map%key(i)
-      prev_key = map%key(i)
-    endif
-  enddo
-  map%n_elements = j
-
-end
-
-subroutine cache_map_shrink(map,thr)
-  use map_module
-  implicit none
-  type (cache_map_type), intent(inout) :: map
-  real(integral_kind)  , intent(in) :: thr
-  integer(cache_map_size_kind)   :: i,j
-
-  j=0
-  do i=1,map%n_elements
-    if (abs(map%value(i)) > thr) then
-      j = j+1
-      map%value(j) = map%value(i)
-      map%key(j) = map%key(i)
-    endif
-  enddo
-  map%n_elements = j
-
-end
-
-subroutine map_unique(map)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  integer(map_size_kind)         :: i
-  integer(map_size_kind)         :: icount
-
-  icount = 0_8
-  !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i)&
-      !$OMP REDUCTION(+:icount)
-  do i=0_8,map%map_size
-    call omp_set_lock(map%map(i)%lock)
-    call cache_map_unique(map%map(i))
-    call omp_unset_lock(map%map(i)%lock)
-    icount = icount + map%map(i)%n_elements
-  enddo
-  !$OMP END PARALLEL DO
-  map%n_elements = icount
-
-end
-
-subroutine map_merge(map)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  integer(map_size_kind)         :: i
-  integer(map_size_kind)         :: icount
-
-  icount = 0_8
-  !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i)&
-      !$OMP REDUCTION(+:icount)
-  do i=0_8,map%map_size
-    call omp_set_lock(map%map(i)%lock)
-    call cache_map_merge(map%map(i))
-    call omp_unset_lock(map%map(i)%lock)
-    icount = icount + map%map(i)%n_elements
-  enddo
-  !$OMP END PARALLEL DO
-  map%n_elements = icount
-
-end
-
-subroutine map_shrink(map,thr)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  real(integral_kind), intent(in) :: thr
-  integer(map_size_kind)         :: i
-  integer(map_size_kind)         :: icount
-
-  icount = 0_8
-  !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i)&
-      !$OMP REDUCTION(+:icount)
-  do i=0_8,map%map_size
-    call omp_set_lock(map%map(i)%lock)
-    call cache_map_shrink(map%map(i),thr)
-    call omp_unset_lock(map%map(i)%lock)
-    icount = icount + map%map(i)%n_elements
-  enddo
-  !$OMP END PARALLEL DO
-  map%n_elements = icount
-
-end
-
-subroutine map_update(map, key, value, sze, thr)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  integer, intent(in)            :: sze
-  integer(key_kind), intent(inout) :: key(sze)
-  real(integral_kind), intent(inout) :: value(sze)
-  real(integral_kind), intent(in) :: thr
-
-  integer                        :: i
-  integer(map_size_kind)         :: idx_cache, idx_cache_new
-  integer(cache_map_size_kind)   :: idx
-  integer                        :: sze2
-  integer(cache_key_kind)        :: cache_key
-  integer(map_size_kind)         :: n_elements_temp
-  type (cache_map_type)          :: local_map
-  logical                        :: map_sorted
-
-  sze2 = sze
-  map_sorted = .True.
-
-  n_elements_temp = 0_8
-  n_elements_temp = n_elements_temp + 1_8
-  do while (sze2>0)
-    i=1
-    do while (i<=sze)
-      if (key(i) /= 0_8) then
-        idx_cache = shiftr(key(i),map_shift)
-        if (omp_test_lock(map%map(idx_cache)%lock)) then
-          local_map%key => map%map(idx_cache)%key
-          local_map%value => map%map(idx_cache)%value
-          local_map%sorted = map%map(idx_cache)%sorted
-          local_map%map_size = map%map(idx_cache)%map_size
-          local_map%n_elements = map%map(idx_cache)%n_elements
-          do
-          !DIR$ FORCEINLINE
-          call search_key_big_interval(key(i),local_map%key, local_map%n_elements, idx, 1, local_map%n_elements)
-          if (idx > 0_8) then
-            local_map%value(idx) = local_map%value(idx) + value(i)
-          else
-            ! Assert that the map has a proper size
-            if (local_map%n_elements == local_map%map_size) then
-              call cache_map_merge(local_map)
-              call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements)
-              call cache_map_shrink(local_map,thr)
-            endif
-            cache_key = int(iand(key(i),map_mask),2)
-            local_map%n_elements = local_map%n_elements + 1
-            local_map%value(local_map%n_elements) = value(i)
-            local_map%key(local_map%n_elements) = cache_key
-            local_map%sorted = .False.
-            n_elements_temp = n_elements_temp + 1_8
-          endif  ! idx > 0
-          key(i) = 0_8
-          i = i+1
-          sze2 = sze2-1
-          if (i>sze) then
-            i=1
-          endif
-          if ( (shiftr(key(i),map_shift) /= idx_cache).or.(key(i)==0_8)) then
-            exit
-          endif
-        enddo
-        map%map(idx_cache)%key => local_map%key
-        map%map(idx_cache)%value => local_map%value
-        map%map(idx_cache)%sorted = local_map%sorted
-        map%map(idx_cache)%n_elements = local_map%n_elements
-        map%map(idx_cache)%map_size = local_map%map_size
-        map_sorted = map_sorted .and. local_map%sorted
-        call omp_unset_lock(map%map(idx_cache)%lock)
-      endif  ! omp_test_lock
-    else
-      i=i+1
-    endif  ! key = 0
-  enddo  ! i
-enddo  ! sze2 > 0
-call omp_set_lock(map%lock)
-map%n_elements = map%n_elements + n_elements_temp
-map%sorted = map%sorted .and. map_sorted
-call omp_unset_lock(map%lock)
-
-end
-
-subroutine map_append(map, key, value, sze)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  integer, intent(in)            :: sze
-  integer(key_kind), intent(inout) :: key(sze)
-  real(integral_kind), intent(inout) :: value(sze)
-
-  integer                        :: i
-  integer(cache_map_size_kind)   :: n_elements
-  integer(map_size_kind)         :: idx_cache
-  integer(cache_key_kind)        :: cache_key
-
-  do i=1,sze
-    idx_cache = shiftr(key(i),map_shift)
-    call omp_set_lock(map%map(idx_cache)%lock)
-    n_elements = map%map(idx_cache)%n_elements + 1
-    ! Assert that the map has a proper size
-    if (n_elements == map%map(idx_cache)%map_size) then
-      call cache_map_reallocate(map%map(idx_cache), n_elements+ shiftr(n_elements,1))
-    endif
-    cache_key = int(iand(key(i),map_mask),2)
-    map%map(idx_cache)%value(n_elements) = value(i)
-    map%map(idx_cache)%key(n_elements) = cache_key
-    map%map(idx_cache)%n_elements = n_elements
-    if (map%map(idx_cache)%sorted.and.n_elements > 1) then
-      map%map(idx_cache)%sorted = (map%map(idx_cache)%key(n_elements-1) <= cache_key)
-      map%sorted = map%sorted .and. map%map(idx_cache)%sorted
-    endif
-    call omp_unset_lock(map%map(idx_cache)%lock)
-  enddo
-  call omp_set_lock(map%lock)
-  map%n_elements = map%n_elements + sze
-  call omp_unset_lock(map%lock)
-
-end
-
-subroutine map_get(map, key, value)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  integer(key_kind), intent(in)  :: key
-  real(integral_kind), intent(out) :: value
-  integer(map_size_kind)         :: idx_cache
-  integer(cache_map_size_kind)   :: idx
-
-  ! index in tha pointers array
-  idx_cache = shiftr(key,map_shift)
-  !DIR$ FORCEINLINE
-  call cache_map_get_interval(map%map(idx_cache), key, value, 1, map%map(idx_cache)%n_elements,idx)
-end
-
-subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx)
-  use map_module
-  implicit none
-  type (cache_map_type), intent(inout) :: map
-  integer(key_kind), intent(in)  :: key
-  integer(cache_map_size_kind), intent(in) :: ibegin, iend
-  real(integral_kind), intent(out) :: value
-  integer(cache_map_size_kind), intent(inout) :: idx
-  double precision, pointer :: v(:)
-  integer :: i
-
-  call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend)
-  if (idx > 0) then
-    value = map%value(idx)
-  else
-    value = 0._integral_kind
-  endif
-!  call search_key_value_big_interval(key, value, map%key, map%value, map%n_elements, idx, ibegin, iend)
-end
-
-
-subroutine map_get_many(map, key, value, sze)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  integer, intent(in)            :: sze
-  integer(key_kind), intent(in)  :: key(sze)
-  real(integral_kind), intent(out) :: value(sze)
-  integer                        :: i
-  integer(map_size_kind)         :: idx_cache
-  integer(cache_map_size_kind)   :: ibegin, iend
-  integer(cache_map_size_kind), allocatable :: idx(:)
-  !DIR$ ATTRIBUTES ALIGN : 64    :: idx
-
-  allocate(idx(sze))
-  do i=1,sze
-    idx_cache = shiftr(key(i),map_shift)
-    iend = map%map(idx_cache)%n_elements
-    !DIR$ FORCEINLINE
-    call search_key_big_interval(key(i),map%map(idx_cache)%key, iend, idx(i), 1, iend)
-  enddo
-  do i=1,sze
-    idx_cache = shiftr(key(i),map_shift)
-    if (idx(i) > 0) then
-      value(i) = map%map(idx_cache)%value(idx(i))
-    else
-      value(i) = 0.
-    endif
-  enddo
-  deallocate(idx)
-end
-
-subroutine map_exists_many(map, key, sze)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  integer, intent(in)            :: sze
-  integer(key_kind), intent(inout) :: key(sze)
-  integer                        :: i
-  integer(map_size_kind)         :: idx_cache, idx_cache_prev
-  integer(cache_map_size_kind)   :: ibegin, iend
-  integer(cache_map_size_kind), allocatable :: idx(:)
-  !DIR$ ATTRIBUTES ALIGN : 64    :: idx
-
-  idx_cache_prev = -1_map_size_kind
-  allocate(idx(sze))
-  do i=1,sze
-    idx_cache = shiftr(key(i),map_shift)
-    iend = map%map(idx_cache)%n_elements
-    if (idx_cache == idx_cache_prev) then
-      if ((idx(i-1) > 0_cache_map_size_kind).and.(idx(i-1) < iend)) then
-        if ((key(i) == key(i-1)+1).and.(map%map(idx_cache)%key(idx(i-1))+1) == key(i)) then
-          idx(i) = idx(i-1)+1
-          cycle
-        endif
-      endif
-    endif
-    !DIR$ FORCEINLINE
-    call search_key_big_interval(key(i),map%map(idx_cache)%key, iend, idx(i), 1, iend)
-    idx_cache_prev = idx_cache
-  enddo
-  do i=1,sze
-    idx_cache = shiftr(key(i),map_shift)
-    if (idx(i) <= 0) then
-      key(i) = 0_key_kind
-    endif
-  enddo
-  deallocate(idx)
-end
-
-subroutine search_key_big(key,X,sze,idx)
-  use map_module
-  implicit none
-  integer(cache_map_size_kind), intent(in) :: sze
-  integer(key_kind)           , intent(in) :: key
-  integer(cache_key_kind)     , intent(in) :: X(sze)
-  integer(cache_map_size_kind), intent(out) :: idx
-
-  call  search_key_big_interval(key,X,sze,idx,1,sze)
-end
-
-
-subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
-  use map_module
-  implicit none
-  integer(cache_map_size_kind), intent(in) :: sze
-  integer(key_kind)           , intent(in) :: key
-  integer(cache_key_kind)     , intent(in) :: X(sze)
-  integer(cache_map_size_kind), intent(in) :: ibegin_in, iend_in
-  integer(cache_map_size_kind), intent(out) :: idx
-
-  integer(cache_map_size_kind)   :: istep, ibegin, iend, i
-  integer(cache_key_kind)        :: cache_key
-
-  if (sze /= 0) then
-    continue
-  else
-    idx = -1
-    return
-  endif
-  cache_key = int(iand(key,map_mask),2)
-  ibegin = min(ibegin_in,sze)
-  iend   = min(iend_in,sze)
-  if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
-
-    istep = shiftr(iend-ibegin,1)
-    idx = ibegin + istep
-    do while (istep > 4)
-      idx = ibegin + istep
-      ! TODO : Cache misses
-      if (cache_key < X(idx)) then
-        iend = idx
-        istep = shiftr(idx-ibegin,1)
-        idx = ibegin + istep
-        if (cache_key < X(idx)) then
-          iend = idx
-          istep = shiftr(idx-ibegin,1)
-          cycle
-        else if (cache_key > X(idx)) then
-          ibegin = idx
-          istep = shiftr(iend-idx,1)
-          cycle
-        else
-          return
-        endif
-      else if (cache_key > X(idx)) then
-        ibegin = idx
-        istep = shiftr(iend-idx,1)
-        idx = idx + istep
-        if (cache_key < X(idx)) then
-          iend = idx
-          istep = shiftr(idx-ibegin,1)
-          cycle
-        else if (cache_key > X(idx)) then
-          ibegin = idx
-          istep = shiftr(iend-idx,1)
-          cycle
-        else
-          return
-        endif
-      else
-        return
-      endif
-    enddo
-    idx = ibegin
-    if (min(iend_in,sze) > ibegin+4) then
-      iend = ibegin+4
-      !DIR$ LOOP COUNT MAX(4)
-      do while (cache_key > X(idx))
-        idx = idx+1
-      end do
-    else
-      !DIR$ LOOP COUNT MAX(4)
-      do while (cache_key > X(idx))
-        idx = idx+1
-        if (idx == iend) then
-          exit
-        endif
-      end do
-    endif
-    if (cache_key /= X(idx)) then
-      idx = 1-idx
-    endif
-    return
-
-  else
-
-    if (cache_key < X(ibegin)) then
-      idx = -ibegin
-      return
-    endif
-    if (cache_key > X(iend)) then
-      idx = -iend
-      return
-    endif
-    if (cache_key == X(ibegin)) then
-      idx = ibegin
-      return
-    endif
-    if (cache_key == X(iend)) then
-      idx = iend
-      return
-    endif
-  endif
-
-end
-
-subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in)
-  use map_module
-  implicit none
-  integer(cache_map_size_kind), intent(in) :: sze
-  integer(key_kind)           , intent(in) :: key
-  real(integral_kind)         , intent(out) :: value
-  integer(cache_key_kind)     , intent(in) :: X(sze)
-  real(integral_kind)         , intent(in) :: Y(sze)
-  integer(cache_map_size_kind), intent(in) :: ibegin_in, iend_in
-  integer(cache_map_size_kind), intent(out) :: idx
-
-  integer(cache_map_size_kind)   :: istep, ibegin, iend, i
-  integer(cache_key_kind)        :: cache_key
-
-  if (sze /= 0) then
-    continue
-  else
-    idx = -1
-    value = 0.d0
-    return
-  endif
-  cache_key = int(iand(key,map_mask),2)
-  ibegin = min(ibegin_in,sze)
-  iend   = min(iend_in,sze)
-  if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
-
-    istep = shiftr(iend+ibegin,1)
-    idx = ibegin + istep
-    do while (istep > 4)
-      idx = ibegin + istep
-      ! TODO : Cache misses
-      if (cache_key < X(idx)) then
-        iend = idx
-        istep = shiftr(idx-ibegin,1)
-        idx = ibegin + istep
-        if (cache_key < X(idx)) then
-          iend = idx
-          istep = shiftr(idx-ibegin,1)
-          cycle
-        else if (cache_key > X(idx)) then
-          ibegin = idx
-          istep = shiftr(iend-idx,1)
-          cycle
-        else
-          value = Y(idx)
-          return
-        endif
-      else if (cache_key > X(idx)) then
-        ibegin = idx
-        istep = shiftr(iend-idx,1)
-        idx = idx + istep
-        if (cache_key < X(idx)) then
-          iend = idx
-          istep = shiftr(idx-ibegin,1)
-          cycle
-        else if (cache_key > X(idx)) then
-          ibegin = idx
-          istep = shiftr(iend-idx,1)
-          cycle
-        else
-          value = Y(idx)
-          return
-        endif
-      else
-        value = Y(idx)
-        return
-      endif
-    enddo
-    idx = ibegin
-    if (min(iend_in,sze) > ibegin+4) then
-      iend = ibegin+4
-      !DIR$ LOOP COUNT MAX(4)
-      do while (cache_key > X(idx))
-        idx = idx+1
-      end do
-    else
-      !DIR$ LOOP COUNT MAX(4)
-      do while (cache_key > X(idx))
-        idx = idx+1
-        if (idx == iend) then
-          exit
-        endif
-      end do
-    endif
-    if (cache_key /= X(idx)) then
-      idx = 1-idx
-      value = 0.d0
-    else
-      value = Y(idx)
-    endif
-    return
-
-  else
-
-    if (cache_key < X(ibegin)) then
-      idx = -ibegin
-      value = 0.d0
-      return
-    endif
-    if (cache_key > X(iend)) then
-      idx = -iend
-      value = 0.d0
-      return
-    endif
-    if (cache_key == X(ibegin)) then
-      idx = ibegin
-      value = Y(idx)
-      return
-    endif
-    if (cache_key == X(iend)) then
-      idx = iend
-      value = Y(idx)
-      return
-    endif
-  endif
-
-end
-
-
-subroutine get_cache_map_n_elements_max(map,n_elements_max)
-  use map_module
-  implicit none
-  ! Returns the size of the largest cache_map
-  type (map_type), intent(in)    :: map
-  integer(cache_map_size_kind), intent(out) :: n_elements_max
-  integer(map_size_kind)         :: i
-  n_elements_max = 0_cache_map_size_kind
-  do i=0_8,map%map_size
-    n_elements_max = max(n_elements_max, map%map(i)%n_elements)
-  enddo
-end
-
-
-
-subroutine get_cache_map(map,map_idx,keys,values,n_elements)
-  use map_module
-  implicit none
-  type (map_type), intent(in)    :: map
-  integer(map_size_kind), intent(in) :: map_idx
-  integer(cache_map_size_kind), intent(inout) :: n_elements
-  integer(key_kind), intent(out) :: keys(n_elements)
-  double precision, intent(out)  :: values(n_elements)
-  integer(cache_map_size_kind)   :: i
-  integer(key_kind)              :: shift
-
-  shift = shiftl(map_idx,map_shift)
-
-  n_elements = map%map(map_idx)%n_elements
-  do i=1,n_elements
-    keys(i)   = map%map(map_idx)%key(i)   + shift
-    values(i) = map%map(map_idx)%value(i)
-  enddo
-
-end
-
diff --git a/src/utils/map_module.mod b/src/utils/map_module.mod
deleted file mode 100644
index a1277b9..0000000
Binary files a/src/utils/map_module.mod and /dev/null differ