diff --git a/src/becke_numerical_grid/grid_becke_vector.irp.f b/src/becke_numerical_grid/grid_becke_vector.irp.f
index 9856bcbd..321d6109 100644
--- a/src/becke_numerical_grid/grid_becke_vector.irp.f
+++ b/src/becke_numerical_grid/grid_becke_vector.irp.f
@@ -5,6 +5,7 @@ BEGIN_PROVIDER [integer, n_points_final_grid]
   END_DOC
   integer                        :: i,j,k,l
   n_points_final_grid = 0
+  print*,'tresh_grid = ',tresh_grid
   do j = 1, nucl_num
     do i = 1, n_points_radial_grid -1
       do k = 1, n_points_integration_angular
@@ -36,6 +37,7 @@ END_PROVIDER
   integer                        :: i,j,k,l,i_count
   double precision               :: r(3)
   i_count = 0
+  print*,'tresh_grid = ',tresh_grid
   do j = 1, nucl_num
     do i = 1, n_points_radial_grid -1
       do k = 1, n_points_integration_angular
diff --git a/src/dft_utils_one_e/sr_coulomb.irp.f b/src/dft_utils_one_e/sr_coulomb.irp.f
index d9d001c4..fea92cc7 100644
--- a/src/dft_utils_one_e/sr_coulomb.irp.f
+++ b/src/dft_utils_one_e/sr_coulomb.irp.f
@@ -35,3 +35,40 @@
  short_range_Hartree = short_range_Hartree * 0.5d0
  print*, 'short_range_Hartree',short_range_Hartree
 END_PROVIDER
+
+
+ BEGIN_PROVIDER [double precision, regular_range_Hartree_operator, (mo_num,mo_num,N_states)]
+&BEGIN_PROVIDER [double precision, regular_range_Hartree, (N_states)]
+ implicit none
+ BEGIN_DOC
+! regular_range_Hartree_operator(i,j) = $\int dr i(r)j(r) \int r' \rho(r') W_{ee}^{sr}$
+!
+! regular_range_Hartree = $1/2  \sum_{i,j} \rho_{ij} \mathtt{regular_range_Hartree_operator}(i,j)$
+!
+!                     = $1/2  \int dr \int r' \rho(r) \rho(r') W_{ee}^{sr}$
+ END_DOC
+ integer :: i,j,k,l,m,n,istate
+ double precision :: get_two_e_integral
+ double precision :: integral, contrib
+ double precision :: integrals_array(mo_num,mo_num)
+ regular_range_Hartree_operator = 0.d0
+ regular_range_Hartree = 0.d0
+ do i = 1, mo_num
+  do j = 1, mo_num
+   if(dabs(one_e_dm_average_mo_for_dft(j,i)).le.1.d-12)cycle
+   call get_mo_two_e_integrals_i1j1(i,j,mo_num,integrals_array,mo_integrals_map)
+   do istate = 1, N_states
+    do k = 1, mo_num
+     do l = 1, mo_num
+      integral = integrals_array(l,k)
+      contrib = one_e_dm_mo_for_dft(i,j,istate) * integral 
+      regular_range_Hartree_operator(l,k,istate) += contrib
+      regular_range_Hartree(istate) += contrib * one_e_dm_mo_for_dft(k,l,istate)
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+ regular_range_Hartree = regular_range_Hartree * 0.5d0
+ print*, 'regular_range_Hartree',regular_range_Hartree
+END_PROVIDER