diff --git a/EZFIO.tar.gz b/EZFIO.tar.gz
new file mode 100644
index 0000000..42dc8d8
Binary files /dev/null and b/EZFIO.tar.gz differ
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..21075c5
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,16 @@
+
+
+
+EZFIO: EZFIO.tar.gz
+	tar -zxvf EZFIO.tar.gz
+	ln -s $$PWD/eplf.config EZFIO/config
+
+EZFIO/lib/libezfio.so: EZFIO eplf.config
+	make -C EZFIO/
+
+bin/ezfio.py: EZFIO/lib/libezfio.so
+	ln -s $$PWD/EZFIO/Python/ezfio.py bin/ezfio.py
+
+bin/eplf: EZFIO/lib/libezfio.so
+	make -C src
+	mv src/eplf bin
diff --git a/README b/README
index cf9815f..e7e91aa 100644
--- a/README
+++ b/README
@@ -5,8 +5,7 @@
 Dependencies:
 -------------
 
-- QCIO library : http://qcio.sourceforge.net
-- resultsFile  : http://resultsFile.sourceforge.net
+- EZFIO library : http://ezfio.sourceforge.net
 - IRPF90 : http://irpf90.sourceforge.net
 - Python > 2.3
 - Any Fortran 90 compiler (gfortran, for example)
diff --git a/src/Makefile b/src/Makefile
index 75e4b3d..aadc8cb 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,11 +1,11 @@
 # MPI-ifort
-IRPF90 = irpf90 -DMPI #-a -d 
+IRPF90 = irpf90 # -DMPI #-a -d 
 FC     = mpif90 -xT -ip -finline
 FCFLAGS= -O3
 
 # Gfortran
 #IRPF90 = irpf90 #-DMPI #-a -d 
-#FC     = gfortran -ffree-line-length-none -static-libgcc 
+#FC     = gfortran -static-libgcc 
 #FCFLAGS= -O3 -ffast-math -L ~/QCIO/lib
 
 # Mono
diff --git a/src/density.irp.f b/src/density.irp.f
index b5dcfd7..742b768 100644
--- a/src/density.irp.f
+++ b/src/density.irp.f
@@ -1,19 +1,101 @@
-BEGIN_PROVIDER [ real, density_p ]
+BEGIN_PROVIDER [ real, density_value_p ]
 
  BEGIN_DOC
-! Value of the density at the current point 
+! Value of the density_value at the current point 
  END_DOC
 
- density_p = 0.
  integer :: i
- do i=1,elec_beta_num
-  density_p = density_p + mo_value_p(i)**2
- enddo
+ density_value_p = density_alpha_value_p + density_beta_value_p
 
- do i=1,elec_alpha_num
-  density_p = density_p + mo_value_p(i)**2
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, density_grad_p, (3) ]
+ implicit none
+ BEGIN_DOC
+! Gradient of the density at the current point 
+ END_DOC
+
+ integer :: i, l
+
+ do l=1,3
+  density_grad_p(l) = density_alpha_grad_p(l) + density_beta_grad_p(l)
  enddo
 
 END_PROVIDER
 
 
+BEGIN_PROVIDER [ real, density_alpha_value_p ]
+
+ BEGIN_DOC
+! Value of the alpha density at the current point 
+ END_DOC
+
+ density_alpha_value_p = 0.
+ integer :: i
+ do i=1,elec_alpha_num
+  density_alpha_value_p = density_alpha_value_p + mo_value_p(i)**2
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ real, density_beta_value_p ]
+
+ BEGIN_DOC
+! Value of the beta density at the current point 
+ END_DOC
+
+ density_beta_value_p = 0.
+ integer :: i
+ do i=1,elec_beta_num
+  density_beta_value_p = density_beta_value_p + mo_value_p(i)**2
+ enddo
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, density_alpha_grad_p, (3) ]
+ implicit none
+ BEGIN_DOC
+! Gradient of the density at the current point 
+ END_DOC
+
+ integer :: i, l
+
+ do l=1,3
+  density_alpha_grad_p(l) = 0.
+ enddo
+
+ do i=1,elec_alpha_num
+  do l=1,3
+   density_alpha_grad_p(l) = density_alpha_grad_p(l) + 2.*mo_grad_p(i,l)*mo_value_p(i)
+  enddo
+ enddo
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, density_beta_grad_p, (3) ]
+ implicit none
+ BEGIN_DOC
+! Gradient of the density at the current point 
+ END_DOC
+
+ integer :: i, l
+
+ do l=1,3
+  density_beta_grad_p(l) = 0.
+ enddo
+
+ do i=1,elec_beta_num
+  do l=1,3
+   density_beta_grad_p(l) = density_beta_grad_p(l) + 2.*mo_grad_p(i,l)*mo_value_p(i)
+  enddo
+ enddo
+
+END_PROVIDER
+
+
+
+
+
diff --git a/src/elf_function.irp.f b/src/elf_function.irp.f
new file mode 100644
index 0000000..e657d9d
--- /dev/null
+++ b/src/elf_function.irp.f
@@ -0,0 +1,67 @@
+BEGIN_PROVIDER [ double precision, kinetic_energy_alpha_p ]
+ implicit none
+ BEGIN_DOC  
+! Kinetic Energy
+ END_DOC
+
+ kinetic_energy_alpha_p = 0.d0
+
+ integer :: i, l
+ do i=1,elec_alpha_num
+  do l=1,3
+   kinetic_energy_alpha_p = kinetic_energy_alpha_p+mo_grad_p(i,l)**2
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, kinetic_energy_beta_p ]
+ implicit none
+ BEGIN_DOC  
+! Kinetic Energy
+ END_DOC
+
+ kinetic_energy_beta_p = 0.d0
+
+ integer :: i,l
+ do i=1,elec_beta_num
+  do l=1,3
+   kinetic_energy_beta_p = kinetic_energy_beta_p+mo_grad_p(i,l)**2
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, kinetic_energy_p ]
+ implicit none
+ BEGIN_DOC  
+ END_DOC
+ kinetic_energy_p = kinetic_energy_alpha_p + kinetic_energy_beta_p
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, elf_value_p ]
+ implicit none
+ BEGIN_DOC  
+! ELF value at point
+ END_DOC
+ double precision, parameter :: Cf = 2.871
+ double precision :: D0, Ds
+ double precision :: modulus_a2, modulus_b2
+
+ integer :: l
+
+ modulus_a2 = 0.d0
+ modulus_b2 = 0.d0
+ do l=1,3
+   modulus_a2 = modulus_a2 + density_alpha_grad_p(l)**2
+   modulus_b2 = modulus_b2 + density_beta_grad_p(l)**2
+ enddo
+
+ Ds = kinetic_energy_p - 0.125d0 *  &
+  ( (modulus_a2/density_alpha_value_p) + &
+    (modulus_b2/density_beta_value_p) )
+ D0 = 2.d0*Cf*density_value_p**(5./3.)
+ elf_value_p = 1.d0/ (1.d0 + (Ds/D0)**2)
+
+END_PROVIDER
diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index 7242946..6fb2253 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -3,9 +3,10 @@ BEGIN_PROVIDER  [ real, eplf_gamma ]
  BEGIN_DOC  
 ! Value of the gaussian for the EPLF
  END_DOC
+ include 'constants.F'
  real            :: eps
  eps = -real(dlog(tiny(1.d0)))
- eplf_gamma = density_p**(2./3.) * eps
+ eplf_gamma = (4./3.*pi*density_value_p)**(2./3.) * eps
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ]
diff --git a/src/mo_point.irp.f b/src/mo_point.irp.f
index f9a3124..b36c1c0 100644
--- a/src/mo_point.irp.f
+++ b/src/mo_point.irp.f
@@ -16,4 +16,40 @@ BEGIN_PROVIDER [ real, mo_value_p, (mo_num) ]
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ real, mo_grad_p, (mo_num,3) ]
+ implicit none
+
+ BEGIN_DOC
+! Gradients of the molecular orbitals
+ END_DOC
+
+ integer :: j, k, l
+
+ do l=1,3
+  do j=1,mo_num
+    mo_grad_p(j,l) = 0.
+    do k=1,ao_num
+      mo_grad_p(j,l) = mo_grad_p(j,l) + mo_coef(k,j)*ao_grad_p(k,l)
+    enddo
+  enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ real, mo_lapl_p, (mo_num) ]
+ implicit none
+
+ BEGIN_DOC
+! Laplacians of the molecular orbitals
+ END_DOC
+
+ integer :: j, k
+
+ do j=1,mo_num
+   mo_lapl_p(j) = 0.
+   do k=1,ao_num
+     mo_lapl_p(j) = mo_lapl_p(j)+mo_coef(k,j)*ao_lapl_p(k)
+   enddo
+ enddo
+
+END_PROVIDER
 
diff --git a/src/mpi.irp.f b/src/mpi.irp.f
index bc0dca5..ed41343 100644
--- a/src/mpi.irp.f
+++ b/src/mpi.irp.f
@@ -16,12 +16,14 @@ subroutine start_mpi
 
 end
 
-BEGIN_PROVIDER [ integer, mpi_rank ]
+BEGIN_PROVIDER [ logical, mpi_master ]
+&BEGIN_PROVIDER [ integer, mpi_rank ]
  implicit none
  BEGIN_DOC  
-! Number of the processor
- END_DOC
+! mpi_rank : Number of the processor
 
+! mpi_master : True if the current processor is the master
+ END_DOC
 
  IRP_IF MPI
    include 'mpif.h'
@@ -38,6 +40,8 @@ BEGIN_PROVIDER [ integer, mpi_rank ]
 
  IRP_ENDIF
 
+ mpi_master = (mpi_rank == 0)
+
 END_PROVIDER
 
 BEGIN_PROVIDER [ integer, mpi_size ]
@@ -65,13 +69,6 @@ BEGIN_PROVIDER [ integer, mpi_size ]
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ logical, mpi_master ]
- implicit none
- BEGIN_DOC  
-! mpi_master : True if the current processor is the master
- END_DOC
- mpi_master = (mpi_rank == 0)
-END_PROVIDER
 
 subroutine barrier()
 IRP_IF MPI
diff --git a/src/test_1d.irp.f b/src/test_1d.irp.f
index 7dfdfec..c617c5c 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_up_up, eplf_up_dn
+  print *,  point(3), eplf_value, eplf_gamma
  enddo
 end
diff --git a/test/c2h.out b/test/c2h.out
index ff7c524..9f19b36 100644
--- a/test/c2h.out
+++ b/test/c2h.out
@@ -1,7 +1,7 @@
  Entering Gaussian System, Link 0=/usr/local/g03/g03
  Initial command:
- /usr/local/g03/l1.exe /tmp/Gau-7519.inp -scrdir=/tmp/
- Entering Link 1 = /usr/local/g03/l1.exe PID=      7523.
+ /usr/local/g03/l1.exe /tmp/Gau-2977.inp -scrdir=/tmp/
+ Entering Link 1 = /usr/local/g03/l1.exe PID=      2981.
   
  Copyright (c) 1988,1990,1992,1993,1995,1998,2003,2004, Gaussian, Inc.
                   All Rights Reserved.
@@ -75,21 +75,21 @@
  
  ******************************************
  Gaussian 03:  AM64L-G03RevC.02 12-Jun-2004
-                28-Mar-2009 
+                29-Oct-2009 
  ******************************************
  %mem=800Mb
  %chk=check
- ---------------------------------------
- #p cc-pvdz ROHF GFPRINT pop=Full 6d 10f
- ---------------------------------------
+ -----------------------------------------------
+ #p cc-pvdz ROHF GFPRINT pop=Full 6d 10f out=wfn
+ -----------------------------------------------
  1/38=1/1;
  2/17=6,18=5,40=1/2;
  3/5=16,8=22,11=2,16=1,24=100,25=1,30=1/1,2,3;
  4/7=6/1;
  5/5=2,32=1,38=5/2;
  6/7=3,28=1/1;
- 99/5=1,9=1/99;
- Leave Link    1 at Sat Mar 28 08:02:08 2009, MaxMem=  104857600 cpu:       0.1
+ 99/5=1,6=100,9=1/99;
+ Leave Link    1 at Thu Oct 29 17:46:24 2009, MaxMem=  104857600 cpu:       0.1
  (Enter /usr/local/g03/l101.exe)
  -----
  titre
@@ -109,7 +109,7 @@
  AtZEff=   0.0000000   0.0000000   0.0000000
  AtQMom=   0.0000000   0.0000000   0.0000000
  AtGFac=   0.0000000   0.0000000   2.7928460
- Leave Link  101 at Sat Mar 28 08:02:09 2009, MaxMem=  104857600 cpu:       0.1
+ Leave Link  101 at Thu Oct 29 17:46:25 2009, MaxMem=  104857600 cpu:       0.1
  (Enter /usr/local/g03/l202.exe)
                           Input orientation:                          
  ---------------------------------------------------------------------
@@ -141,7 +141,7 @@
     3          1             0        0.000000    0.000000   -1.531957
  ---------------------------------------------------------------------
  Rotational constants (GHZ):      0.0000000     44.9353403     44.9353403
- Leave Link  202 at Sat Mar 28 08:02:10 2009, MaxMem=  104857600 cpu:       0.1
+ Leave Link  202 at Thu Oct 29 17:46:26 2009, MaxMem=  104857600 cpu:       0.1
  (Enter /usr/local/g03/l301.exe)
  Standard basis: CC-pVDZ (6D, 10F)
  AO basis set:
@@ -219,17 +219,17 @@
  ScaDFX=  1.000000  1.000000  1.000000  1.000000
  IRadAn=      0 IRanWt=     -1 IRanGd=            0 ICorTp=0
  NAtoms=    3 NActive=    3 NUniq=    3 SFac= 1.00D+00 NAtFMM=   60 Big=F
- Leave Link  301 at Sat Mar 28 08:02:11 2009, MaxMem=  104857600 cpu:       0.0
+ Leave Link  301 at Thu Oct 29 17:46:27 2009, MaxMem=  104857600 cpu:       0.0
  (Enter /usr/local/g03/l302.exe)
  NPDir=0 NMtPBC=     1 NCelOv=     1 NCel=       1 NClECP=     1 NCelD=      1
          NCelK=      1 NCelE2=     1 NClLst=     1 CellRange=     0.0.
  One-electron integrals computed using PRISM.
  NBasis=    35 RedAO= T  NBF=    19     2     7     7
  NBsUse=    35 1.00D-06 NBFU=    19     2     7     7
- Leave Link  302 at Sat Mar 28 08:02:13 2009, MaxMem=  104857600 cpu:       0.1
+ Leave Link  302 at Thu Oct 29 17:46:28 2009, MaxMem=  104857600 cpu:       0.1
  (Enter /usr/local/g03/l303.exe)
  DipDrv:  MaxL=1.
- Leave Link  303 at Sat Mar 28 08:02:14 2009, MaxMem=  104857600 cpu:       0.0
+ Leave Link  303 at Thu Oct 29 17:46:29 2009, MaxMem=  104857600 cpu:       0.0
  (Enter /usr/local/g03/l401.exe)
  Harris functional with IExCor=  205 diagonalized for initial guess.
  ExpMin= 1.22D-01 ExpMax= 6.67D+03 ExpMxC= 2.28D+02 IAcc=1 IRadAn=         1 AccDes= 1.00D-06
@@ -243,7 +243,7 @@
                  (DLTA) (PI) (PI) (SG) (SG) (PI) (PI) (SG) (SG)
                  (SG)
  The electronic state of the initial guess is 2-SG.
- Leave Link  401 at Sat Mar 28 08:02:16 2009, MaxMem=  104857600 cpu:       0.0
+ Leave Link  401 at Thu Oct 29 17:46:30 2009, MaxMem=  104857600 cpu:       0.0
  (Enter /usr/local/g03/l502.exe)
  Warning!  Cutoffs for single-point calculations used.
  Restricted open shell SCF:
@@ -339,7 +339,7 @@
  KE= 7.596465821717D+01 PE=-2.179438569834D+02 EE= 4.551191533841D+01
  Annihilation of the first spin contaminant:
  S**2 before annihilation     0.7500,   after     0.7500
- Leave Link  502 at Sat Mar 28 08:02:17 2009, MaxMem=  104857600 cpu:       0.2
+ Leave Link  502 at Thu Oct 29 17:46:31 2009, MaxMem=  104857600 cpu:       0.2
  (Enter /usr/local/g03/l601.exe)
  Copying SCF densities to generalized density rwf, ISCF=0 IROHF=1.
 
@@ -1245,20 +1245,20 @@
  ---------------------------------------------------------------------------------
 
  No NMR shielding tensors so no spin-rotation constants.
- Leave Link  601 at Sat Mar 28 08:02:19 2009, MaxMem=  104857600 cpu:       0.3
+ Leave Link  601 at Thu Oct 29 17:46:33 2009, MaxMem=  104857600 cpu:       0.3
  (Enter /usr/local/g03/l9999.exe)
- 1\1\GINC-LPQSV11\SP\ROHF\CC-pVDZ\C2H1(2)\SCEMAMA\28-Mar-2009\0\\#P CC-
- PVDZ ROHF GFPRINT POP=FULL 6D 10F\\titre\\0,2\C,0,0.,0.,0.59801\C,0,0.
- ,0.,-0.59801\H,0,0.,0.,1.65962\\Version=AM64L-G03RevC.02\State=2-SG\HF
- =-76.1419687\RMSD=8.827e-05\Dipole=0.,0.,0.2930001\PG=C*V [C*(H1C1C1)]
- \\@
+ 1\1\GINC-LPQSV11\SP\ROHF\CC-pVDZ\C2H1(2)\SCEMAMA\29-Oct-2009\0\\#P CC-
+ PVDZ ROHF GFPRINT POP=FULL 6D 10F OUT=WFN\\titre\\0,2\C,0,0.,0.,0.5980
+ 1\C,0,0.,0.,-0.59801\H,0,0.,0.,1.65962\\Version=AM64L-G03RevC.02\State
+ =2-SG\HF=-76.1419687\RMSD=8.827e-05\Dipole=0.,0.,0.2930001\PG=C*V [C*(
+ H1C1C1)]\\@
+
+ Writing a WFN file to c2h.wfn.
 
 
- IT WAS AN ACT OF DESPARATION. FOR SIX YEARS I HAD STRUGGLED WITH THE
- BLACKBODY THEORY. I KNEW THE PROBLEM WAS FUNDAMENTAL, AND I KNEW
- THE ANSWER. I HAD TO FIND A THEORETICAL EXPLANATION AT ANY COST,
- EXCEPT FOR THE INVIOLABLITY OF THE TWO LAWS OF THERMODYNAMICS.
-   -- MAX PLANCK, 1931
- Job cpu time:  0 days  0 hours  0 minutes  9.0 seconds.
- File lengths (MBytes):  RWF=     13 Int=      0 D2E=      0 Chk=     11 Scr=      1
- Normal termination of Gaussian 03 at Sat Mar 28 08:02:20 2009.
+ IN THE WOODS WE RETURN TO REASON AND FAITH.
+
+                                 -- EMERSON
+ Job cpu time:  0 days  0 hours  0 minutes  8.8 seconds.
+ File lengths (MBytes):  RWF=     13 Int=      0 D2E=      0 Chk=     10 Scr=      1
+ Normal termination of Gaussian 03 at Thu Oct 29 17:46:34 2009.