diff --git a/eplf.config b/eplf.config
index 1de93f5..ca7f5b5 100644
--- a/eplf.config
+++ b/eplf.config
@@ -37,4 +37,29 @@ grid
  step_size           real (3)
  origin              real (3)
  opposite            real (3)
+ num_x               integer           = grid_point_num(1)
+ num_y               integer           = grid_point_num(2)
+ num_z               integer           = grid_point_num(3)
+
+grid_data
+ eplf                real (grid_num_x,grid_num_y,grid_num_z)
+ eplf_grad           real (grid_num_x,grid_num_y,grid_num_z,4)
+ eplf_lapl           real (grid_num_x,grid_num_y,grid_num_z)
+ elf                 real (grid_num_x,grid_num_y,grid_num_z)
+ elf_grad            real (grid_num_x,grid_num_y,grid_num_z,4)
+ elf_lapl            real (grid_num_x,grid_num_y,grid_num_z)
+ density             real (grid_num_x,grid_num_y,grid_num_z)
+ density_grad        real (grid_num_x,grid_num_y,grid_num_z,4)
+ density_lapl        real (grid_num_x,grid_num_y,grid_num_z)
+
+compute
+ eplf                logical
+ eplf_grad           logical
+ eplf_lapl           logical
+ elf                 logical
+ elf_grad            logical
+ elf_lapl            logical
+ density             logical
+ density_grad        logical
+ density_lapl        logical
 
diff --git a/src/Makefile b/src/Makefile
index aadc8cb..df8a1e4 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,5 +1,5 @@
 # MPI-ifort
-IRPF90 = irpf90 # -DMPI #-a -d 
+IRPF90 = irpf90 -DMPI #-a -d 
 FC     = mpif90 -xT -ip -finline
 FCFLAGS= -O3
 
@@ -17,6 +17,9 @@ SRC=
 OBJ=
 LIB=../EZFIO/lib/libezfio.a
 
+eplf: main
+	mv main eplf
+
 include irpf90.make
 
 irpf90.make: $(wildcard *.irp.f)
diff --git a/src/compute.irp.f b/src/compute.irp.f
new file mode 100644
index 0000000..f2e6d86
--- /dev/null
+++ b/src/compute.irp.f
@@ -0,0 +1,29 @@
+BEGIN_SHELL [ /usr/bin/python ]
+
+to_compute = [\
+ "eplf",
+ "eplf_grad",
+ "eplf_lapl",
+ "elf",
+ "elf_grad",
+ "elf_lapl",
+ "density",
+ "density_grad",
+ "density_lapl",
+]
+
+template = """
+BEGIN_PROVIDER [ logical, comp_$X ]
+ implicit none
+ BEGIN_DOC  
+! If true, $X
+ END_DOC
+ comp_$X = .False.
+ call get_compute_$X(comp_$X)
+END_PROVIDER
+"""
+
+for t in to_compute:
+  print template.replace("$X",t)
+
+END_SHELL
diff --git a/src/debug_eplf.irp.f b/src/debug_eplf.irp.f
index 1ab649b..a02f679 100644
--- a/src/debug_eplf.irp.f
+++ b/src/debug_eplf.irp.f
@@ -46,8 +46,8 @@ program debug
  print *, 'EPLF integral N    :', ao_eplf_integral_numeric(i,j,eplf_gamma,point)
  
  print *, ''
- print *, 'EPLF grid Npoints  :', grid_eplf_x_num, grid_eplf_y_num, grid_eplf_z_num
- print *, 'EPLF grid step     :', grid_eplf_step(:)
- print *, 'EPLF grid origin   :', grid_eplf_origin(:)
+ print *, 'EPLF grid Npoints  :', grid_x_num, grid_y_num, grid_z_num
+ print *, 'EPLF grid step     :', grid_step(:)
+ print *, 'EPLF grid origin   :', grid_origin(:)
 
 end
diff --git a/src/eplf.irp.f b/src/eplf.irp.f
deleted file mode 100644
index c764287..0000000
--- a/src/eplf.irp.f
+++ /dev/null
@@ -1,6 +0,0 @@
-program eplf
- provide mpi_rank
- call write_grid_eplf()
- call finish()
-end
-
diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index 6fb2253..a7088d4 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -111,7 +111,7 @@ END_PROVIDER
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ real, eplf_value ]
+BEGIN_PROVIDER [ real, eplf_value_p ]
  implicit none
  BEGIN_DOC  
 ! Value of the EPLF at the current point.
@@ -128,9 +128,9 @@ BEGIN_PROVIDER [ real, eplf_value ]
    ab = -(dlog(ab)/eplf_gamma)
    aa = dsqrt(aa)
    ab = dsqrt(ab)
-   eplf_value = (aa-ab)/(aa+ab+eps)
+   eplf_value_p = (aa-ab)/(aa+ab+eps)
  else
-   eplf_value = 0.d0
+   eplf_value_p = 0.d0
  endif
 
 END_PROVIDER
diff --git a/src/eplf_grid.irp.f b/src/eplf_grid.irp.f
deleted file mode 100644
index db17b3c..0000000
--- a/src/eplf_grid.irp.f
+++ /dev/null
@@ -1,163 +0,0 @@
-BEGIN_PROVIDER [ character*(100), grid_cube_filename ]
- BEGIN_DOC  
-! Name of the file containing the parameters of the grid
- END_DOC
- grid_cube_filename = 'eplf_grid.cube'
-END_PROVIDER
-
- BEGIN_PROVIDER [ integer, grid_eplf_x_num ]
-&BEGIN_PROVIDER [ integer, grid_eplf_y_num ]
-&BEGIN_PROVIDER [ integer, grid_eplf_z_num ]
-&BEGIN_PROVIDER [ real   , grid_eplf_step  , (3) ]
-&BEGIN_PROVIDER [ real   , grid_eplf_origin, (3) ]
-
- real, parameter :: UNDEFINED=123456789.123456789
-
- BEGIN_DOC  
-! Number of grid points in x, y, z directions
- END_DOC
-
- integer :: Npoints(3)
- real :: step_size(3)
- real :: origin(3)
- real :: opposite(3)
-
- Npoints  (:)  = 80
- origin   (:)  = UNDEFINED
- opposite (:)  = UNDEFINED
- step_size(:)  = UNDEFINED
-
- call get_grid_point_num(Npoints)
- call get_grid_origin(origin)
- call get_grid_opposite(opposite)
- call get_grid_step_size(step_size)
- 
- if (origin(1) == UNDEFINED) then
-   integer :: i,l
-   do l=1,3
-    origin(l) = nucl_coord(1,l)
-    do i=2,nucl_num
-     origin(l) = min(origin(l),nucl_coord(i,l))
-    enddo
-    origin(l) = origin(l) - 4.
-   enddo
- endif
-
- if (opposite(1) == UNDEFINED) then
-   do l=1,3
-    opposite(l) = nucl_coord(1,l)
-    do i=2,nucl_num
-     opposite(l) = max(opposite(l),nucl_coord(i,l))
-    enddo
-    opposite(l) = opposite(l) + 4.
-   enddo
- endif
-
- if (step_size(1) == UNDEFINED) then
-   do l=1,3
-     step_size(l) = (opposite(l) - origin(l))/float(Npoints(l))
-   enddo
- endif
-
- do l=1,3
-   grid_eplf_origin(l) = origin(l)
-   grid_eplf_step(l)   = step_size(l)
- enddo
- grid_eplf_x_num = Npoints(1)
- grid_eplf_y_num = Npoints(2)
- grid_eplf_z_num = Npoints(3)
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ real, grid_eplf, (grid_eplf_x_num,grid_eplf_y_num,grid_eplf_z_num) ]
- implicit none
- BEGIN_DOC
-! EPLF on a grid
- END_DOC
-
- IRP_IF MPI
-   include 'mpif.h'
- IRP_ENDIF
-
- integer :: ix, iy, iz
- integer :: ibegin, iend
-
- do iz=1,grid_eplf_z_num 
-  do iy=1,grid_eplf_y_num 
-   do ix=1,grid_eplf_x_num 
-    grid_eplf(ix,iy,iz) = 0.
-   enddo
-  enddo
- enddo
-
- integer :: icount
- icount = mpi_size
- do iz=1,grid_eplf_z_num
-  if (mpi_master) then
-    print *,  int(100*dble(iz)/dble(grid_eplf_z_num)), '%'
-  endif
-  point(3) = grid_eplf_origin(3)+(iz-1)*grid_eplf_step(3)
-  do iy=1,grid_eplf_y_num 
-   point(2) = grid_eplf_origin(2)+(iy-1)*grid_eplf_step(2)
-   do ix=1,grid_eplf_x_num 
-    icount = icount-1
-    if (icount == mpi_rank) then
-     point(1) = grid_eplf_origin(1)+(ix-1)*grid_eplf_step(1)
-     TOUCH point
-     grid_eplf(ix,iy,iz) = eplf_value
-    endif
-    if (icount == 0) then
-      icount = mpi_size
-    endif
-   enddo
-  enddo
- enddo
-
- IRP_IF MPI
-   integer :: dim, ierr
-   do iz=1,grid_eplf_z_num 
-    real   :: buffer(grid_eplf_x_num*grid_eplf_y_num)
-    icount = 0
-    do iy=1,grid_eplf_y_num 
-     do ix=1,grid_eplf_x_num 
-      buffer(icount+ix) = grid_eplf(ix,iy,iz)
-     enddo
-     icount = icount + grid_eplf_x_num
-    enddo
-    dim = grid_eplf_x_num * grid_eplf_y_num
-    call MPI_REDUCE(buffer,grid_eplf(1,1,iz),dim,mpi_real, &
-       mpi_sum,0,MPI_COMM_WORLD,ierr)
-   enddo
- IRP_ENDIF
-
-END_PROVIDER
-
-subroutine write_grid_eplf
- implicit none
- integer :: i
- integer :: l
- integer :: ix, iy, iz
- if (.not.mpi_master) then
-  return
- endif
- open(unit=99,file=grid_cube_filename,status='UNKNOWN',action='WRITE')
- write (99,*) 'Cube File'
- write (99,*) 'Analytical EPLF grid'
- write (99,10) nucl_num,(grid_eplf_origin(i), i=1,3)
- write (99,10) grid_eplf_x_num, grid_eplf_step(1), 0., 0.
- write (99,10) grid_eplf_y_num, 0., grid_eplf_step(2), 0.
- write (99,10) grid_eplf_z_num, 0., 0., grid_eplf_step(3)
- do i=1,nucl_num
-  write (99,11) int(nucl_charge(i)), nucl_charge(i), (nucl_coord(i,l),l=1,3)
- enddo
- do ix = 1, grid_eplf_x_num
-  do iy = 1, grid_eplf_y_num
-    write (99,20) (grid_eplf(ix,iy,iz), iz=1, grid_eplf_z_num)
-  enddo
- enddo
- 10   format (2X,I3,3(2X,F10.6))
- 11   format (2X,I3,4(2X,F10.6))
- 20   format (6(E13.5))
- close(99)
-end    
-
diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f
index d0aa953..ea44b6f 100644
--- a/src/ezfio_interface.irp.f
+++ b/src/ezfio_interface.irp.f
@@ -21,6 +21,15 @@ data = [ \
 ("grid_step_size"                  , "real"          , "(3)"                   ),
 ("grid_origin"                     , "real"          , "(3)"                   ),
 ("grid_opposite"                   , "real"          , "(3)"                   ),
+("compute_eplf"                    , "logical"       , ""                      ),
+("compute_eplf_grad"               , "logical"       , ""                      ),
+("compute_eplf_lapl"               , "logical"       , ""                      ),
+("compute_elf"                     , "logical"       , ""                      ),
+("compute_elf_grad"                , "logical"       , ""                      ),
+("compute_elf_lapl"                , "logical"       , ""                      ),
+("compute_density"                 , "logical"       , ""                      ),
+("compute_density_grad"            , "logical"       , ""                      ),
+("compute_density_lapl"            , "logical"       , ""                      ),
 ]
 
 data_no_set = [\
diff --git a/src/finish.irp.f b/src/finish.irp.f
index 93e5a40..c74f602 100644
--- a/src/finish.irp.f
+++ b/src/finish.irp.f
@@ -5,14 +5,12 @@ subroutine abrt (here,message)
   IRP_ENDIF
   character*(*) :: here
   character*(*) :: message
-  if (mpi_master) then
-   print *, ''
-   print *, '-------------------------'
-   print *, 'Error in '//trim(here)//':'
-   print *, '-------------------------'
-   print *, trim(message)//'.'
-   print *, '-------------------------'
-  endif
+  print *, ''
+  print *, '-------------------------'
+  print *, 'Error in '//trim(here)//':'
+  print *, '-------------------------'
+  print *, trim(message)//'.'
+  print *, '-------------------------'
 
   IRP_IF MPI
    integer :: ierr
diff --git a/src/gradients.irp.f b/src/gradients.irp.f
new file mode 100644
index 0000000..978fa97
--- /dev/null
+++ b/src/gradients.irp.f
@@ -0,0 +1,41 @@
+BEGIN_SHELL [ /usr/bin/python ]
+
+values = [\
+ "eplf",
+ "elf",
+ ]
+
+template = """
+ BEGIN_PROVIDER [ real, $X_grad_p, (3) ]
+&BEGIN_PROVIDER [ real, $X_lapl_p ]
+ implicit none
+ BEGIN_DOC  
+! Gradient and Laplacian of the EPLF at the current point.
+ END_DOC
+
+ real, parameter :: Delta=0.001
+ integer :: l
+
+ $X_lapl_p = -6.*$X_value_p
+
+ do l=1,3
+  point(l) = point(l)+Delta
+  TOUCH point
+  $X_grad_p(l) = $X_value_p
+  $X_lapl_p    = $X_lapl_p + $X_value_p
+  point(l) = point(l)-Delta-Delta
+  TOUCH point
+  $X_grad_p(l) = 0.5*($X_grad_p(l) - $X_value_p)/Delta
+  $X_lapl_p    = $X_lapl_p + $X_value_p
+  point(l) = point(l) + Delta
+ enddo
+ $X_lapl_p = $X_lapl_p/Delta**2
+ TOUCH point
+
+END_PROVIDER
+"""
+
+for value in values:
+  print template.replace("$X",value)
+
+END_SHELL
diff --git a/src/grid.irp.f b/src/grid.irp.f
new file mode 100644
index 0000000..fb4f8cd
--- /dev/null
+++ b/src/grid.irp.f
@@ -0,0 +1,284 @@
+ BEGIN_PROVIDER [ integer, grid_x_num ]
+&BEGIN_PROVIDER [ integer, grid_y_num ]
+&BEGIN_PROVIDER [ integer, grid_z_num ]
+&BEGIN_PROVIDER [ real   , grid_step  , (3) ]
+&BEGIN_PROVIDER [ real   , grid_origin, (3) ]
+
+ real, parameter :: UNDEFINED=123456789.123456789
+
+ BEGIN_DOC  
+! Number of grid points in x, y, z directions
+ END_DOC
+
+ integer :: Npoints(3)
+ real :: step_size(3)
+ real :: origin(3)
+ real :: opposite(3)
+
+ Npoints  (:)  = 80
+ origin   (:)  = UNDEFINED
+ opposite (:)  = UNDEFINED
+ step_size(:)  = UNDEFINED
+
+ call get_grid_point_num(Npoints)
+ call get_grid_origin(origin)
+ call get_grid_opposite(opposite)
+ call get_grid_step_size(step_size)
+ 
+ if (origin(1) == UNDEFINED) then
+   integer :: i,l
+   do l=1,3
+    origin(l) = nucl_coord(1,l)
+    do i=2,nucl_num
+     origin(l) = min(origin(l),nucl_coord(i,l))
+    enddo
+    origin(l) = origin(l) - 4.
+   enddo
+ endif
+
+ if (opposite(1) == UNDEFINED) then
+   do l=1,3
+    opposite(l) = nucl_coord(1,l)
+    do i=2,nucl_num
+     opposite(l) = max(opposite(l),nucl_coord(i,l))
+    enddo
+    opposite(l) = opposite(l) + 4.
+   enddo
+ endif
+
+ if (step_size(1) == UNDEFINED) then
+   do l=1,3
+     step_size(l) = (opposite(l) - origin(l))/float(Npoints(l))
+   enddo
+ endif
+
+ do l=1,3
+   grid_origin(l) = origin(l)
+   grid_step(l)   = step_size(l)
+ enddo
+ grid_x_num = Npoints(1)
+ grid_y_num = Npoints(2)
+ grid_z_num = Npoints(3)
+
+ call get_grid_point_num(Npoints)
+ call get_grid_origin(origin)
+ call get_grid_opposite(opposite)
+ call get_grid_step_size(step_size)
+
+END_PROVIDER
+
+!subroutine write_grid_eplf
+! implicit none
+! integer :: i
+! integer :: l
+! integer :: ix, iy, iz
+! if (.not.mpi_master) then
+!  return
+! endif
+! open(unit=99,file=grid_cube_filename,status='UNKNOWN',action='WRITE')
+! write (99,*) 'Cube File'
+! write (99,*) 'Analytical EPLF grid'
+! write (99,10) nucl_num,(grid_eplf_origin(i), i=1,3)
+! write (99,10) grid_eplf_x_num, grid_eplf_step(1), 0., 0.
+! write (99,10) grid_eplf_y_num, 0., grid_eplf_step(2), 0.
+! write (99,10) grid_eplf_z_num, 0., 0., grid_eplf_step(3)
+! do i=1,nucl_num
+!  write (99,11) int(nucl_charge(i)), nucl_charge(i), (nucl_coord(i,l),l=1,3)
+! enddo
+! do ix = 1, grid_eplf_x_num
+!  do iy = 1, grid_eplf_y_num
+!    write (99,20) (grid_eplf(ix,iy,iz), iz=1, grid_eplf_z_num)
+!  enddo
+! enddo
+! 10   format (2X,I3,3(2X,F10.6))
+! 11   format (2X,I3,4(2X,F10.6))
+! 20   format (6(E13.5))
+! close(99)
+!end    
+
+BEGIN_SHELL [ /usr/bin/python ]
+
+grids = [ \
+ "eplf",
+ "elf",
+ "density", 
+]
+
+template = """
+BEGIN_PROVIDER [ real, grid_$X, (grid_x_num,grid_y_num,grid_z_num) ]
+ implicit none
+ BEGIN_DOC
+! $X on a grid
+ END_DOC
+
+ IRP_IF MPI
+   include 'mpif.h'
+ IRP_ENDIF
+
+ integer :: ix, iy, iz
+ integer :: ibegin, iend
+
+ do iz=1,grid_z_num 
+  do iy=1,grid_y_num 
+   do ix=1,grid_x_num 
+    grid_$X(ix,iy,iz) = 0.
+   enddo
+  enddo
+ enddo
+
+ integer :: icount
+ icount = mpi_size
+ do iz=1,grid_z_num
+  if (mpi_master) then
+    print *,  int(100*dble(iz)/dble(grid_z_num)), '%'
+  endif
+  point(3) = grid_origin(3)+(iz-1)*grid_step(3)
+  do iy=1,grid_y_num 
+   point(2) = grid_origin(2)+(iy-1)*grid_step(2)
+   do ix=1,grid_x_num 
+    icount = icount-1
+    if (icount == mpi_rank) then
+     point(1) = grid_origin(1)+(ix-1)*grid_step(1)
+     TOUCH point
+     grid_$X(ix,iy,iz) = $X_value_p
+    endif
+    if (icount == 0) then
+      icount = mpi_size
+    endif
+   enddo
+  enddo
+ enddo
+
+ IRP_IF MPI
+   integer :: dim, ierr
+   do iz=1,grid_z_num 
+    real   :: buffer(grid_x_num*grid_y_num)
+    icount = 0
+    do iy=1,grid_y_num 
+     do ix=1,grid_x_num 
+      buffer(icount+ix) = grid_$X(ix,iy,iz)
+     enddo
+     icount = icount + grid_x_num
+    enddo
+    dim = grid_x_num * grid_y_num
+    call MPI_REDUCE(buffer,grid_$X(1,1,iz),dim,mpi_real, &
+       mpi_sum,0,MPI_COMM_WORLD,ierr)
+   enddo
+ IRP_ENDIF
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ real, grid_$X_grad, (grid_x_num,grid_y_num,grid_z_num,4) ]
+&BEGIN_PROVIDER [ real, grid_$X_lapl, (grid_x_num,grid_y_num,grid_z_num) ]
+ BEGIN_DOC
+! Laplacian of $X on a grid
+ END_DOC
+ implicit none
+ BEGIN_DOC
+! Gradient and lapacian of $X on a grid. 4th dimension of the grad is its norm.
+ END_DOC
+
+ IRP_IF MPI
+   include 'mpif.h'
+ IRP_ENDIF
+
+ integer :: ix, iy, iz, it
+ integer :: ibegin, iend
+
+ do iz=1,grid_z_num 
+  do iy=1,grid_y_num 
+   do ix=1,grid_x_num 
+    do it=1,4
+     grid_$X_grad(ix,iy,iz,it) = 0.
+    enddo
+    grid_$X_lapl(ix,iy,iz) = 0.
+   enddo
+  enddo
+ enddo
+
+ integer :: icount
+ icount = mpi_size
+ do iz=1,grid_z_num
+  if (mpi_master) then
+    print *,  int(100*dble(iz)/dble(grid_z_num)), '%'
+  endif
+  point(3) = grid_origin(3)+(iz-1)*grid_step(3)
+  do iy=1,grid_y_num 
+   point(2) = grid_origin(2)+(iy-1)*grid_step(2)
+   do ix=1,grid_x_num 
+    icount = icount-1
+    if (icount == mpi_rank) then
+     point(1) = grid_origin(1)+(ix-1)*grid_step(1)
+     TOUCH point
+     do it=1,3
+      grid_$X_grad(ix,iy,iz,it) = $X_grad_p(it)
+     enddo
+     grid_$X_grad(ix,iy,iz,4) = $X_grad_p(1)**2 + $X_grad_p(2)**2 + $X_grad_p(3)**2
+     grid_$X_lapl(ix,iy,iz) = $X_lapl_p
+     grid_$X_grad(ix,iy,iz,4) = sqrt(grid_$X_grad(ix,iy,iz,4))
+    endif
+    if (icount == 0) then
+      icount = mpi_size
+    endif
+   enddo
+  enddo
+ enddo
+
+ IRP_IF MPI
+  integer :: dim, ierr
+  do it=1,4 
+   do iz=1,grid_z_num 
+    real   :: buffer(grid_x_num*grid_y_num)
+    icount = 0
+    do iy=1,grid_y_num 
+     do ix=1,grid_x_num 
+      buffer(icount+ix) = grid_$X_grad(ix,iy,iz,it)
+     enddo
+     icount = icount + grid_x_num
+    enddo
+    dim = grid_x_num * grid_y_num
+    call MPI_REDUCE(buffer,grid_$X_grad(1,1,iz,it),dim,mpi_real, &
+       mpi_sum,0,MPI_COMM_WORLD,ierr)
+    icount = 0
+    do iy=1,grid_y_num 
+     do ix=1,grid_x_num 
+      buffer(icount+ix) = grid_$X_lapl(ix,iy,iz)
+     enddo
+     icount = icount + grid_x_num
+    enddo
+    dim = grid_x_num * grid_y_num
+    call MPI_REDUCE(buffer,grid_$X_lapl(1,1,iz),dim,mpi_real, &
+       mpi_sum,0,MPI_COMM_WORLD,ierr)
+   enddo
+  enddo
+ IRP_ENDIF
+
+END_PROVIDER
+
+subroutine set_grid_data_$X(buffer)
+ real   :: buffer(grid_x_num,grid_y_num,grid_z_num)
+ if (mpi_master) then
+   call ezfio_set_grid_data_$X(buffer)
+ endif
+end
+
+subroutine set_grid_data_$X_grad(buffer)
+ real   :: buffer(grid_x_num,grid_y_num,grid_z_num,4)
+ if (mpi_master) then
+   call ezfio_set_grid_data_$X_grad(buffer)
+ endif
+end
+
+subroutine set_grid_data_$X_lapl(buffer)
+ real   :: buffer(grid_x_num,grid_y_num,grid_z_num)
+ if (mpi_master) then
+   call ezfio_set_grid_data_$X_lapl(buffer)
+ endif
+end
+"""
+
+for grid in grids: 
+  print template.replace("$X",grid)
+
+END_SHELL
diff --git a/src/main.irp.f b/src/main.irp.f
new file mode 100644
index 0000000..b38f416
--- /dev/null
+++ b/src/main.irp.f
@@ -0,0 +1,39 @@
+program eplf
+ implicit none
+ provide mpi_rank
+
+ if (comp_eplf) then
+   call set_grid_data_eplf(grid_eplf)
+   FREE grid_eplf
+ endif
+ if (comp_eplf_grad.or.comp_eplf_lapl) then
+   call set_grid_data_eplf_grad(grid_eplf_grad)
+   call set_grid_data_eplf_lapl(grid_eplf_lapl)
+   FREE grid_eplf_grad
+   FREE grid_eplf_lapl
+ endif
+
+ if (comp_elf) then
+   call set_grid_data_elf(grid_elf)
+   FREE grid_elf
+ endif
+ if (comp_elf_grad.or.comp_elf_lapl) then
+   call set_grid_data_elf_grad(grid_elf_grad)
+   call set_grid_data_elf_lapl(grid_elf_lapl)
+   FREE grid_elf_grad
+   FREE grid_elf_lapl
+ endif
+
+ if (comp_density) then
+   call set_grid_data_density(grid_density)
+   FREE grid_density
+ endif
+ if (comp_density_grad.or.comp_density_lapl) then
+   call set_grid_data_density_grad(grid_density_grad)
+   call set_grid_data_density_lapl(grid_density_lapl)
+   FREE grid_density_grad
+   FREE grid_density_lapl
+ endif
+ call finish()
+end
+
diff --git a/src/mpi.irp.f b/src/mpi.irp.f
index ed41343..63459ac 100644
--- a/src/mpi.irp.f
+++ b/src/mpi.irp.f
@@ -1,60 +1,26 @@
-subroutine start_mpi
- implicit none
- integer :: ierr
- integer, save :: started = 0
-
- IRP_IF MPI
-   include 'mpif.h'
-   if (started == 0) then
-     call MPI_INIT(ierr)
-     if (ierr /= MPI_SUCCESS) then
-       call abrt(irp_here,"Unable to initialize MPI")
-     endif
-   endif
-   started = 1
- IRP_ENDIF
-
-end
-
 BEGIN_PROVIDER [ logical, mpi_master ]
 &BEGIN_PROVIDER [ integer, mpi_rank ]
+&BEGIN_PROVIDER [ integer, mpi_size ]
  implicit none
  BEGIN_DOC  
 ! mpi_rank : Number of the processor
 
+! mpi_size : Number of processors
+
 ! mpi_master : True if the current processor is the master
  END_DOC
 
  IRP_IF MPI
    include 'mpif.h'
    integer :: ierr
-   call start_mpi
+   call MPI_INIT(ierr)
+   if (ierr /= MPI_SUCCESS) then
+     call abrt(irp_here,"Unable to initialize MPI")
+   endif
    call MPI_COMM_RANK(MPI_COMM_WORLD, mpi_rank, ierr)
    if (ierr /= MPI_SUCCESS) then
-     call abrt(irp_here,"Unable to get MPI")
+     call abrt(irp_here,"Unable to get MPI rank")
    endif
-
- IRP_ELSE
-
-   mpi_rank = 0
-
- IRP_ENDIF
-
- mpi_master = (mpi_rank == 0)
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer, mpi_size ]
- implicit none
- BEGIN_DOC  
-! Number of processors
- END_DOC
-
- IRP_IF MPI
-   include 'mpif.h'
-
-   integer :: ierr
-   call start_mpi
    call MPI_COMM_SIZE(MPI_COMM_WORLD, mpi_size, ierr)
    if (ierr /= MPI_SUCCESS) then
      call abrt(irp_here,"Unable to get MPI size")
@@ -62,10 +28,13 @@ BEGIN_PROVIDER [ integer, mpi_size ]
 
  IRP_ELSE
 
+   mpi_rank = 0
    mpi_size = 1
 
  IRP_ENDIF
 
+ mpi_master = (mpi_rank == 0)
+
 END_PROVIDER
 
 
diff --git a/src/test_1d.irp.f b/src/test_1d.irp.f
index c617c5c..f10c38a 100644
--- a/src/test_1d.irp.f
+++ b/src/test_1d.irp.f
@@ -13,6 +13,6 @@ subroutine run
  do i=- 60,40
   point(3) = real(i)/10.
   TOUCH point
-  print *,  point(3), eplf_value, eplf_gamma
+  print *,  point(3), eplf_value_p, eplf_gamma
  enddo
 end