diff --git a/external/irpf90 b/external/irpf90
index 43160c60..4ab1b175 160000
--- a/external/irpf90
+++ b/external/irpf90
@@ -1 +1 @@
-Subproject commit 43160c60d88d9f61fb97cc0b35477c8eb0df862b
+Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6
diff --git a/plugins/local/extra_basis_int/NEED b/plugins/local/extra_basis_int/NEED
index 0419f87f..76ae72cd 100644
--- a/plugins/local/extra_basis_int/NEED
+++ b/plugins/local/extra_basis_int/NEED
@@ -2,3 +2,4 @@ ao_extra_basis
 ao_one_e_ints
 ao_two_e_ints
 determinants
+becke_numerical_grid
diff --git a/plugins/local/extra_basis_int/test_extra_basis.irp.f b/plugins/local/extra_basis_int/test_extra_basis.irp.f
index b394c47a..71663767 100644
--- a/plugins/local/extra_basis_int/test_extra_basis.irp.f
+++ b/plugins/local/extra_basis_int/test_extra_basis.irp.f
@@ -26,12 +26,21 @@ program test_extra_basis
    print*,j,i,pot_vne_A_extra_basis(j,i)
   enddo
  enddo
- print*,'Density matrix from system A'
+ print*,'Density matrix from system A in the AO basis'
  do i = 1, ao_num
-  do j = 1, ao_num
-   print*,j,i,one_e_dm_ao(j,i)
-  enddo
+!  do j = 1, ao_num
+   write(*,'(100(F16.10,X))')one_e_dm_ao(:,i)
+!  enddo
  enddo
+
+ print*,'Density matrix from system A in the MO basis'
+ do i = 1, mo_num
+!  do j = 1, mo_num
+   write(*,'(100(F16.10,X))')one_e_dm_mo(:,i)
+!  enddo
+ enddo
+
+
  print*,'Density matrix from system B'
  output=trim(ezfio_filename)//'.one_e_dm_b'
  i_unit_output = getUnitAndOpen(output,'w')
diff --git a/src/ao_extra_basis/NEED b/src/ao_extra_basis/NEED
index c4e8c3cf..1680d011 100644
--- a/src/ao_extra_basis/NEED
+++ b/src/ao_extra_basis/NEED
@@ -1,3 +1,4 @@
 extra_nuclei
 basis
 ao_basis
+becke_numerical_grid
diff --git a/src/ao_extra_basis/qp_copy_extra_basis b/src/ao_extra_basis/qp_copy_extra_basis
index 6d0e17b8..9d08392b 100755
--- a/src/ao_extra_basis/qp_copy_extra_basis
+++ b/src/ao_extra_basis/qp_copy_extra_basis
@@ -10,7 +10,8 @@ dir=${QP}
 EZFIO_extra=${1%/}
 EZFIO_extra=${EZFIO_extra%.xyz}
 EZFIO_target=${2%/}
-
+qp set_file ${EZFIO_extra}
+qp run write_extra_grid
 
 
 echo "********** SCRIPT TO COPY DATA FROM EZFIO TO ANOTHER *********"
@@ -65,3 +66,24 @@ i=data_one_e_dm_tot_ao.gz
 newfile=ao_extra_one_e_dm.gz
 cp ${EZFIO_extra}/aux_quantities/$i ${EZFIO_target}/ao_extra_basis/$newfile
 
+echo "COPYING ALL DATA FROM "$EZFIO_extra" extra grid to "${EZFIO_target}"/becke_numerical_grid/"
+## NPOINTS
+i=n_points_extra_final_grid
+newfile=n_points_extra_final_grid
+cp ${EZFIO_extra}/becke_numerical_grid/$i ${EZFIO_target}/becke_numerical_grid/$newfile
+
+## WEIGHTS
+i=final_weight_at_r_vector_extra.gz
+newfile=final_weight_at_r_vector_extra.gz
+cp ${EZFIO_extra}/becke_numerical_grid/$i ${EZFIO_target}/becke_numerical_grid/$newfile
+
+## GRID POINTS 
+i=final_grid_points_extra.gz
+newfile=final_grid_points_extra.gz
+cp ${EZFIO_extra}/becke_numerical_grid/$i ${EZFIO_target}/becke_numerical_grid/$newfile
+
+## INPUT/OUTPUT 
+i=io_extra_grid
+newfile=io_extra_grid
+cp ${EZFIO_extra}/becke_numerical_grid/$i ${EZFIO_target}/becke_numerical_grid/$newfile
+
diff --git a/src/ao_extra_basis/write_extra_grid.irp.f b/src/ao_extra_basis/write_extra_grid.irp.f
new file mode 100644
index 00000000..ffeee5b7
--- /dev/null
+++ b/src/ao_extra_basis/write_extra_grid.irp.f
@@ -0,0 +1,13 @@
+program write_extra_grid_in_ezfio
+ implicit none
+ io_extra_grid = "Write"
+ touch io_extra_grid
+ call routine
+
+end
+
+subroutine routine
+ implicit none
+ provide final_grid_points_extra
+ call ezfio_set_becke_numerical_grid_io_extra_grid("Read")
+end
diff --git a/src/becke_numerical_grid/EZFIO.cfg b/src/becke_numerical_grid/EZFIO.cfg
index e660fd6d..74973e34 100644
--- a/src/becke_numerical_grid/EZFIO.cfg
+++ b/src/becke_numerical_grid/EZFIO.cfg
@@ -38,6 +38,18 @@ type: integer
 doc: Total number of extra_grid points 
 interface: ezfio
 
+[final_grid_points_extra]
+type: double precision
+doc: Grid points on the extra grid
+size: (3,becke_numerical_grid.n_points_extra_final_grid)
+interface: ezfio
+
+[final_weight_at_r_vector_extra]
+type: double precision
+doc: Weights of the points on the extra grid
+size: (becke_numerical_grid.n_points_extra_final_grid)
+interface: ezfio
+
 [extra_grid_type_sgn]
 type: integer
 doc: Type of extra_grid used for the Becke's numerical extra_grid. Can be, by increasing accuracy: [ 0 | 1 | 2 | 3 ]
@@ -80,3 +92,9 @@ doc: method used to sample the radial space. Possible choices are [KNOWLES | GIL
 interface: ezfio,provider,ocaml
 default: KNOWLES
 
+[io_extra_grid]
+type: Disk_access
+doc: Read/Write the extra grid from/to disk [ Write | Read | None ]
+interface: ezfio,provider,ocaml
+default: None
+
diff --git a/src/becke_numerical_grid/extra_grid_vector.irp.f b/src/becke_numerical_grid/extra_grid_vector.irp.f
index e054e22c..bb0296a5 100644
--- a/src/becke_numerical_grid/extra_grid_vector.irp.f
+++ b/src/becke_numerical_grid/extra_grid_vector.irp.f
@@ -10,22 +10,30 @@ BEGIN_PROVIDER [integer, n_points_extra_final_grid]
   implicit none
   integer :: i, j, k, l
 
-  n_points_extra_final_grid = 0
-
-  do j = 1, nucl_num
-    do i = 1, n_points_extra_radial_grid -1
-      do k = 1, n_points_extra_integration_angular
-        if(dabs(final_weight_at_r_extra(k,i,j)) < thresh_extra_grid) then
-          cycle
-        endif
-        n_points_extra_final_grid += 1
+  if (read_extra_grid)then
+   print*,'Reading number of extra grid points '
+   call ezfio_get_becke_numerical_grid_n_points_extra_final_grid(i)
+   n_points_extra_final_grid = i
+  else
+    print*,'Computing number of extra grid points and weights based on a sperical grid'
+    n_points_extra_final_grid = 0
+    do j = 1, nucl_num
+      do i = 1, n_points_extra_radial_grid -1
+        do k = 1, n_points_extra_integration_angular
+          if(dabs(final_weight_at_r_extra(k,i,j)) < thresh_extra_grid) then
+            cycle
+          endif
+          n_points_extra_final_grid += 1
+        enddo
       enddo
     enddo
-  enddo
-
+   print*, ' n max point               = ', n_points_extra_integration_angular*(n_points_extra_radial_grid*nucl_num - 1)
+  endif
   print*, ' n_points_extra_final_grid = ', n_points_extra_final_grid
-  print*, ' n max point               = ', n_points_extra_integration_angular*(n_points_extra_radial_grid*nucl_num - 1)
-  call ezfio_set_becke_numerical_grid_n_points_extra_final_grid(n_points_extra_final_grid)
+  if(write_extra_grid)then
+   print*,'Writing number of extra grid points and weights based on a sperical grid'
+   call ezfio_set_becke_numerical_grid_n_points_extra_final_grid(n_points_extra_final_grid)
+  endif
 
 END_PROVIDER
 
@@ -54,26 +62,41 @@ END_PROVIDER
   call wall_time(wall0)
   print *, ' Providing extra_final_grid_points ...'
 
-  i_count = 0
-  do j = 1, nucl_num
-    do i = 1, n_points_extra_radial_grid -1
-      do k = 1, n_points_extra_integration_angular
-        if(dabs(final_weight_at_r_extra(k,i,j)) < thresh_extra_grid)then
-          cycle
-        endif
-        i_count += 1
-        final_grid_points_extra(1,i_count) = grid_points_extra_per_atom(1,k,i,j)
-        final_grid_points_extra(2,i_count) = grid_points_extra_per_atom(2,k,i,j)
-        final_grid_points_extra(3,i_count) = grid_points_extra_per_atom(3,k,i,j)
-        final_weight_at_r_vector_extra(i_count) = final_weight_at_r_extra(k,i,j)
-        index_final_points_extra(1,i_count) = k
-        index_final_points_extra(2,i_count) = i
-        index_final_points_extra(3,i_count) = j
-        index_final_points_extra_reverse(k,i,j) = i_count
-      enddo
-    enddo
-  enddo
+ 
+  if (read_extra_grid)then
+   print*,'Reading extra grid points and weights'
+   call ezfio_get_becke_numerical_grid_final_grid_points_extra(final_grid_points_extra)
+   call ezfio_get_becke_numerical_grid_final_weight_at_r_vector_extra(final_weight_at_r_vector_extra)
+   index_final_points_extra = -10000000
+   index_final_points_extra_reverse = -10000
+  else
+   print*,'Computing extra grid points and weights based on a sperical grid'
+   i_count = 0
+   do j = 1, nucl_num
+     do i = 1, n_points_extra_radial_grid -1
+       do k = 1, n_points_extra_integration_angular
+         if(dabs(final_weight_at_r_extra(k,i,j)) < thresh_extra_grid)then
+           cycle
+         endif
+         i_count += 1
+         final_grid_points_extra(1,i_count) = grid_points_extra_per_atom(1,k,i,j)
+         final_grid_points_extra(2,i_count) = grid_points_extra_per_atom(2,k,i,j)
+         final_grid_points_extra(3,i_count) = grid_points_extra_per_atom(3,k,i,j)
+         final_weight_at_r_vector_extra(i_count) = final_weight_at_r_extra(k,i,j)
+         index_final_points_extra(1,i_count) = k
+         index_final_points_extra(2,i_count) = i
+         index_final_points_extra(3,i_count) = j
+         index_final_points_extra_reverse(k,i,j) = i_count
+       enddo
+     enddo
+   enddo
+  endif
 
+  if(write_extra_grid)then
+   print*,'Writing extra grid points and weights based on a sperical grid'
+   call ezfio_set_becke_numerical_grid_final_grid_points_extra(final_grid_points_extra)
+   call ezfio_set_becke_numerical_grid_final_weight_at_r_vector_extra(final_weight_at_r_vector_extra)
+  endif
   call wall_time(wall1)
   print *, ' wall time for extra_final_grid_points,', wall1 - wall0
   call print_memory_usage()
diff --git a/src/becke_numerical_grid/grid_becke_vector.irp.f b/src/becke_numerical_grid/grid_becke_vector.irp.f
index f696b023..e9d4e3cf 100644
--- a/src/becke_numerical_grid/grid_becke_vector.irp.f
+++ b/src/becke_numerical_grid/grid_becke_vector.irp.f
@@ -22,7 +22,8 @@ BEGIN_PROVIDER [integer, n_points_final_grid]
 
   print*,' n_points_final_grid = ', n_points_final_grid
   print*,' n max point         = ', n_points_integration_angular*(n_points_radial_grid*nucl_num - 1)
-  call ezfio_set_becke_numerical_grid_n_points_final_grid(n_points_final_grid)
+ ! no reason to write in the EZFIO file the number of grid points ?
+!  call ezfio_set_becke_numerical_grid_n_points_final_grid(n_points_final_grid)
 
 END_PROVIDER