diff --git a/external/ezfio b/external/ezfio
index ed1df9f3..ccee52d0 160000
--- a/external/ezfio
+++ b/external/ezfio
@@ -1 +1 @@
-Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c
+Subproject commit ccee52d00c2cde1d628b0d34f4a247143747bf36
diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg
index 1c66e758..42123373 100644
--- a/src/basis/EZFIO.cfg
+++ b/src/basis/EZFIO.cfg
@@ -15,15 +15,10 @@ interface: ezfio, provider
 
 [shell_normalization_factor]
 type: double precision
-doc: Number of primitives per |AO|
+doc: Normalization factor applied to the whole shell, ex $1/\sqrt{ <d_{z^2}|d_{z^2}>}$
 size: (basis.shell_num)
 interface: ezfio
 
-[prim_num]
-type: integer
-doc: Total number of primitives
-interface: ezfio, provider
-
 [shell_ang_mom]
 type: integer
 doc: Angular momentum of each shell
@@ -48,13 +43,24 @@ doc: Index of the nucleus on which the shell is centered
 size: (basis.shell_num)
 interface: ezfio, provider
 
-[shell_prim_coef]
+[prim_normalization_factor]
+type: double precision
+doc: Normalization factor applied to each primitive
+size: (basis.prim_num)
+interface: ezfio
+
+[prim_num]
+type: integer
+doc: Total number of primitives
+interface: ezfio, provider
+
+[prim_coef]
 type: double precision
 doc: Primitive coefficients
 size: (basis.prim_num)
 interface: ezfio, provider
 
-[shell_prim_expo]
+[prim_expo]
 type: double precision
 doc: Exponents in the shell
 size: (basis.prim_num)
diff --git a/src/basis/basis.irp.f b/src/basis/basis.irp.f
index e25e6a46..a247a90e 100644
--- a/src/basis/basis.irp.f
+++ b/src/basis/basis.irp.f
@@ -15,7 +15,6 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ]
       call ezfio_get_basis_shell_normalization_factor(shell_normalization_factor)
     else
 
-
       double precision               :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
       integer                        :: l, powA(3), nz
       integer                        :: i,j,k
@@ -25,23 +24,22 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ]
       C_A(3) = 0.d0
 
       do i=1,shell_num
+
         powA(1) = shell_ang_mom(i)
         powA(2) = 0
         powA(3) = 0
 
-        ! Normalization of the contracted basis functions
         norm = 0.d0
-        do j=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1
-          do k=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1
-            call overlap_gaussian_xyz(C_A,C_A,shell_prim_expo(j),shell_prim_expo(k),&
-                powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
-            norm = norm+c*shell_prim_coef(j)*shell_prim_coef(k)
+        do k=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1
+          do j=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1
+            call overlap_gaussian_xyz(C_A,C_A,prim_expo(j),prim_expo(k), &
+              powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
+            norm = norm+c*prim_coef(j)*prim_coef(k)
           enddo
         enddo
-        shell_normalization_factor(i) = dsqrt(norm)
+        shell_normalization_factor(i) = 1.d0/dsqrt(norm)
       enddo
 
-
     endif
   endif
   IRP_IF MPI_DEBUG
@@ -60,3 +58,62 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ]
   call write_time(6)
 
 END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Number of primitives per |AO|
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  if (mpi_master) then
+    if (size(prim_normalization_factor) == 0) return
+
+    call ezfio_has_basis_prim_normalization_factor(has)
+    if (has) then
+      write(6,'(A)') '.. >>>>> [ IO READ: prim_normalization_factor ] <<<<< ..'
+      call ezfio_get_basis_prim_normalization_factor(prim_normalization_factor)
+    else
+
+      double precision               :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
+      integer                        :: l, powA(3), nz
+      integer                        :: i,j,k
+      nz=100
+      C_A(1) = 0.d0
+      C_A(2) = 0.d0
+      C_A(3) = 0.d0
+
+      do i=1,shell_num
+
+        powA(1) = shell_ang_mom(i)
+        powA(2) = 0
+        powA(3) = 0
+
+        do k=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1
+            call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), &
+              powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
+            prim_normalization_factor(k) = 1.d0/dsqrt(norm)
+        enddo
+      enddo
+
+
+    endif
+  endif
+  IRP_IF MPI_DEBUG
+  print *,  irp_here, mpi_rank
+  call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+  include 'mpif.h'
+  integer                        :: ierr
+  call MPI_BCAST( prim_normalization_factor, (prim_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+  if (ierr /= MPI_SUCCESS) then
+    stop 'Unable to read prim_normalization_factor with MPI'
+  endif
+  IRP_ENDIF
+
+  call write_time(6)
+
+END_PROVIDER