diff --git a/plugins/local/basis_correction/pbe_on_top.irp.f b/plugins/local/basis_correction/pbe_on_top.irp.f
index be3a23d7..96d1ffdc 100644
--- a/plugins/local/basis_correction/pbe_on_top.irp.f
+++ b/plugins/local/basis_correction/pbe_on_top.irp.f
@@ -180,4 +180,29 @@
  enddo
  END_PROVIDER
 
-
+ BEGIN_PROVIDER [ double precision, extrapolated_on_top, (n_points_final_grid,N_states)]
+&BEGIN_PROVIDER [ double precision, average_extrapolated_on_top, (N_states)]
+ implicit none
+ double precision :: weight,on_top, on_top_extrap, mu_correction_of_on_top,mu
+ integer :: istate, ipoint
+ extrapolated_on_top = 0.d0
+ average_extrapolated_on_top = 0.d0
+ do istate = 1, N_states
+  do ipoint = 1, n_points_final_grid
+   weight = final_weight_at_r_vector(ipoint)
+   mu    = mu_of_r_prov(ipoint,istate)
+   if(mu_of_r_potential == "cas_full")then
+    ! You take the on-top of the CAS wave function which is computed with mu(r) 
+    on_top = on_top_cas_mu_r(ipoint,istate)
+   else
+    ! You take the on-top of the CAS wave function computed separately
+    on_top = total_cas_on_top_density(ipoint,istate)
+   endif
+!  We take the extrapolated on-top pair density 
+   on_top_extrap = mu_correction_of_on_top(mu,on_top)
+   extrapolated_on_top(ipoint,istate) = on_top_extrap
+   average_extrapolated_on_top(istate) += on_top_extrap * weight
+  enddo
+ enddo
+ 
+ END_PROVIDER 
diff --git a/plugins/local/basis_correction/print_n2_related_stuff.irp.f b/plugins/local/basis_correction/print_n2_related_stuff.irp.f
new file mode 100644
index 00000000..0967d16e
--- /dev/null
+++ b/plugins/local/basis_correction/print_n2_related_stuff.irp.f
@@ -0,0 +1,17 @@
+program print_n2_stuffs
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+  read_wf = .True.
+  touch read_wf
+  no_core_density = .True.
+  touch no_core_density
+  call routine
+end
+
+subroutine routine
+ implicit none
+ print*,'average_extrapolated_on_top = ',average_extrapolated_on_top
+ print*,'average_on_top              = ',average_on_top
+end
diff --git a/src/cas_based_on_top/cas_based_on_top.irp.f b/src/cas_based_on_top/cas_based_on_top.irp.f
index 89516336..7d9090ce 100644
--- a/src/cas_based_on_top/cas_based_on_top.irp.f
+++ b/src/cas_based_on_top/cas_based_on_top.irp.f
@@ -16,4 +16,5 @@ end
 subroutine routine
  implicit none
  provide total_cas_on_top_density
+ provide average_on_top
 end
diff --git a/src/cas_based_on_top/on_top_cas_prov.irp.f b/src/cas_based_on_top/on_top_cas_prov.irp.f
index 18a7b7d2..dd93ed40 100644
--- a/src/cas_based_on_top/on_top_cas_prov.irp.f
+++ b/src/cas_based_on_top/on_top_cas_prov.irp.f
@@ -74,3 +74,17 @@
 
  END_PROVIDER 
 
+
+ BEGIN_PROVIDER [ double precision, average_on_top, (n_states)]
+ implicit none
+ integer :: i_point,istate
+ double precision :: wall_0,wall_1,core_inact_dm,pure_act_on_top_of_r,weight
+ average_on_top = 0.d0
+ do istate = 1, N_states
+  do i_point = 1, n_points_final_grid
+   weight = final_weight_at_r_vector(i_point)
+   average_on_top(istate) += total_cas_on_top_density(i_point,istate) * weight
+  enddo
+ enddo
+ print*,'Average on top pair density = ',average_on_top
+ END_PROVIDER