diff --git a/src/dft_utils_in_r/dm_in_r.irp.f b/src/dft_utils_in_r/dm_in_r.irp.f
index 6a19fed0..18eb5403 100644
--- a/src/dft_utils_in_r/dm_in_r.irp.f
+++ b/src/dft_utils_in_r/dm_in_r.irp.f
@@ -109,6 +109,101 @@ end
    grad_dm_b *= 2.d0
  end
 
+
+
+ subroutine density_and_grad_lapl_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, grad_dm_a, grad_dm_b, lapl_dm_a, lapl_dm_b, aos_array, grad_aos_array, lapl_aos_array)
+ implicit none
+ BEGIN_DOC
+! input:
+!
+! * r(1) ==> r(1) = x, r(2) = y, r(3) = z
+!
+! output:
+!
+! * dm_a = alpha density evaluated at r
+! * dm_b = beta  density evaluated at r
+! * aos_array(i) = ao(i) evaluated at r
+! * grad_dm_a(1) = X gradient of the alpha density evaluated in r
+! * grad_dm_a(1) = X gradient of the beta  density evaluated in r
+! * grad_aos_array(1) = X gradient of the aos(i) evaluated at r
+!
+ END_DOC
+ double precision, intent(in)  :: r(3)
+ double precision, intent(out) :: dm_a(N_states),dm_b(N_states)
+ double precision, intent(out) :: grad_dm_a(3,N_states),grad_dm_b(3,N_states)
+ double precision, intent(out) :: lapl_dm_a(3,N_states),lapl_dm_b(3,N_states)
+ double precision, intent(out) :: grad_aos_array(3,ao_num)
+ double precision, intent(out) :: lapl_aos_array(3,ao_num)
+ integer :: i,j,istate
+ double precision  :: aos_array(ao_num),aos_array_bis(ao_num),u_dot_v
+ double precision  :: aos_grad_array(ao_num,3), aos_grad_array_bis(ao_num,3)
+ double precision  :: aos_lapl_array(ao_num,3)
+
+ call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,grad_aos_array,lapl_aos_array)
+ do i = 1, ao_num
+  do j = 1, 3
+   aos_grad_array(i,j) =  grad_aos_array(j,i)
+   aos_lapl_array(i,j) =  lapl_aos_array(j,i)
+  enddo
+ enddo
+
+ do istate = 1, N_states
+   ! alpha density
+   ! aos_array_bis = \rho_ao * aos_array
+   call dsymv('U',ao_num,1.d0,one_e_dm_alpha_ao_for_dft(1,1,istate),size(one_e_dm_alpha_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
+   dm_a(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
+
+   ! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i)
+   grad_dm_a(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
+   grad_dm_a(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
+   grad_dm_a(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
+
+   ! lapl_dm(1) = \sum_i aos_lapl_array(i,1) * aos_array_bis(i)
+   lapl_dm_a(1,istate) = 2.d0 * u_dot_v(aos_lapl_array(1,1),aos_array_bis,ao_num)
+   lapl_dm_a(2,istate) = 2.d0 * u_dot_v(aos_lapl_array(1,2),aos_array_bis,ao_num)
+   lapl_dm_a(3,istate) = 2.d0 * u_dot_v(aos_lapl_array(1,3),aos_array_bis,ao_num)
+
+   ! aos_grad_array_bis(1) = \rho_ao * aos_grad_array(1)
+   call dsymv('U',ao_num,1.d0,one_e_dm_alpha_ao_for_dft(1,1,istate),size(one_e_dm_alpha_ao_for_dft,1),aos_grad_array(1,1),1,0.d0,aos_grad_array_bis(1,1),1)
+   call dsymv('U',ao_num,1.d0,one_e_dm_alpha_ao_for_dft(1,1,istate),size(one_e_dm_alpha_ao_for_dft,1),aos_grad_array(1,2),1,0.d0,aos_grad_array_bis(1,2),1)
+   call dsymv('U',ao_num,1.d0,one_e_dm_alpha_ao_for_dft(1,1,istate),size(one_e_dm_alpha_ao_for_dft,1),aos_grad_array(1,3),1,0.d0,aos_grad_array_bis(1,3),1)
+   ! lapl_dm(1) += \sum_i aos_grad_array(i,1) * aos_grad_array_bis(i)
+   lapl_dm_a(1,istate) += 2.d0 * u_dot_v(aos_grad_array(1,1),aos_grad_array_bis,ao_num)
+   lapl_dm_a(2,istate) += 2.d0 * u_dot_v(aos_grad_array(1,2),aos_grad_array_bis,ao_num)
+   lapl_dm_a(3,istate) += 2.d0 * u_dot_v(aos_grad_array(1,3),aos_grad_array_bis,ao_num)
+ 
+
+   ! beta density
+   call dsymv('U',ao_num,1.d0,one_e_dm_beta_ao_for_dft(1,1,istate),size(one_e_dm_beta_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
+   dm_b(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
+
+   ! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i)
+   grad_dm_b(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
+   grad_dm_b(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
+   grad_dm_b(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
+
+   ! lapl_dm(1) = \sum_i aos_lapl_array(i,1) * aos_array_bis(i)
+   lapl_dm_b(1,istate) = 2.d0 * u_dot_v(aos_lapl_array(1,1),aos_array_bis,ao_num)
+   lapl_dm_b(2,istate) = 2.d0 * u_dot_v(aos_lapl_array(1,2),aos_array_bis,ao_num)
+   lapl_dm_b(3,istate) = 2.d0 * u_dot_v(aos_lapl_array(1,3),aos_array_bis,ao_num)
+
+   ! aos_grad_array_bis(1) = \rho_ao * aos_grad_array(1)
+   call dsymv('U',ao_num,1.d0,one_e_dm_alpha_ao_for_dft(1,1,istate),size(one_e_dm_alpha_ao_for_dft,1),aos_grad_array(1,1),1,0.d0,aos_grad_array_bis(1,1),1)
+   call dsymv('U',ao_num,1.d0,one_e_dm_alpha_ao_for_dft(1,1,istate),size(one_e_dm_alpha_ao_for_dft,1),aos_grad_array(1,2),1,0.d0,aos_grad_array_bis(1,2),1)
+   call dsymv('U',ao_num,1.d0,one_e_dm_alpha_ao_for_dft(1,1,istate),size(one_e_dm_alpha_ao_for_dft,1),aos_grad_array(1,3),1,0.d0,aos_grad_array_bis(1,3),1)
+   ! lapl_dm(1) += \sum_i aos_grad_array(i,1) * aos_grad_array_bis(i)
+   lapl_dm_b(1,istate) += 2.d0 * u_dot_v(aos_grad_array(1,1),aos_grad_array_bis,ao_num)
+   lapl_dm_b(2,istate) += 2.d0 * u_dot_v(aos_grad_array(1,2),aos_grad_array_bis,ao_num)
+   lapl_dm_b(3,istate) += 2.d0 * u_dot_v(aos_grad_array(1,3),aos_grad_array_bis,ao_num)
+ enddo
+   grad_dm_a *= 2.d0
+   grad_dm_b *= 2.d0
+
+ end
+
+
+
+
 subroutine dm_dft_alpha_beta_no_core_at_r(r,dm_a,dm_b)
  implicit none
  BEGIN_DOC
@@ -257,12 +352,13 @@ END_PROVIDER
 &BEGIN_PROVIDER [double precision, one_e_dm_and_grad_beta_in_r,  (4,n_points_final_grid,N_states) ]
 &BEGIN_PROVIDER [double precision, one_e_grad_2_dm_alpha_at_r, (n_points_final_grid,N_states) ]
 &BEGIN_PROVIDER [double precision, one_e_grad_2_dm_beta_at_r, (n_points_final_grid,N_states) ]
+&BEGIN_PROVIDER [double precision, one_e_grad_dm_squared_at_r, (3,n_points_final_grid,N_states) ]
  BEGIN_DOC
 ! one_e_dm_and_grad_alpha_in_r(1,i,i_state) = d\dx n_alpha(r_i,istate)
 ! one_e_dm_and_grad_alpha_in_r(2,i,i_state) = d\dy n_alpha(r_i,istate)
 ! one_e_dm_and_grad_alpha_in_r(3,i,i_state) = d\dz n_alpha(r_i,istate)
 ! one_e_dm_and_grad_alpha_in_r(4,i,i_state) = n_alpha(r_i,istate)
-! one_e_grad_2_dm_alpha_at_r(i,istate)      = d\dx n_alpha(r_i,istate)^2 + d\dy n_alpha(r_i,istate)^2 + d\dz n_alpha(r_i,istate)^2
+! one_e_grad_2_dm_alpha_at_r(i,istate)      = (d\dx n_alpha(r_i,istate))^2 + (d\dy n_alpha(r_i,istate))^2 + (d\dz n_alpha(r_i,istate))^2
 ! where r_i is the ith point of the grid and istate is the state number
  END_DOC
  implicit none
@@ -291,6 +387,9 @@ END_PROVIDER
    one_e_dm_and_grad_beta_in_r(3,i,istate)  =  dm_b_grad(3,istate)
    one_e_dm_and_grad_beta_in_r(4,i,istate)  =  dm_b(istate)
    one_e_grad_2_dm_beta_at_r(i,istate) = dm_b_grad(1,istate) * dm_b_grad(1,istate) + dm_b_grad(2,istate) * dm_b_grad(2,istate) + dm_b_grad(3,istate) * dm_b_grad(3,istate)
+   one_e_grad_dm_squared_at_r(1,i,istate) = 2.D0 * (dm_a_grad(1,istate) + dm_b_grad(1,istate) ) * (one_e_dm_and_grad_alpha_in_r(4,i,istate) + one_e_dm_and_grad_beta_in_r(4,i,istate))
+   one_e_grad_dm_squared_at_r(2,i,istate) = 2.D0 * (dm_a_grad(2,istate) + dm_b_grad(2,istate) ) * (one_e_dm_and_grad_alpha_in_r(4,i,istate) + one_e_dm_and_grad_beta_in_r(4,i,istate))
+   one_e_grad_dm_squared_at_r(3,i,istate) = 2.D0 * (dm_a_grad(3,istate) + dm_b_grad(3,istate) ) * (one_e_dm_and_grad_alpha_in_r(4,i,istate) + one_e_dm_and_grad_beta_in_r(4,i,istate))
   enddo
  enddo
 
diff --git a/src/dft_utils_one_e/ec_lyp.irp.f b/src/dft_utils_one_e/ec_lyp.irp.f
new file mode 100644
index 00000000..0e3daa65
--- /dev/null
+++ b/src/dft_utils_one_e/ec_lyp.irp.f
@@ -0,0 +1,121 @@
+subroutine give_all_stuffs_in_r_for_lyp_88(r,rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_2)
+ implicit none
+ double precision, intent(in) :: r(3)
+ double precision, intent(out) :: rho_a(N_states),rho_b(N_states),grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_2(N_states),rho(N_states)
+ double precision :: grad_rho_a(3,N_states),grad_rho_b(3,N_states),grad_rho_a_b(N_states)
+ double precision :: grad_aos_array(3,ao_num),aos_array(ao_num)
+
+ call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,rho_a,rho_b, grad_rho_a, grad_rho_b, aos_array, grad_aos_array)
+ integer :: i,istate
+ rho = rho_a + rho_b
+ grad_rho_a_2 = 0.d0
+ grad_rho_b_2 = 0.d0
+ grad_rho_a_b = 0.d0
+ do istate = 1, N_states
+  do i = 1, 3 
+   grad_rho_a_2(istate) += grad_rho_a(i,istate) * grad_rho_a(i,istate)
+   grad_rho_b_2(istate) += grad_rho_b(i,istate) * grad_rho_b(i,istate)
+   grad_rho_a_b(istate) += grad_rho_a(i,istate) * grad_rho_b(i,istate)
+  enddo
+ enddo
+ grad_rho_2 = grad_rho_a_2 + grad_rho_b_2 + 2.d0 * grad_rho_a_b
+
+end
+
+
+double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_2)
+
+ implicit none
+
+ BEGIN_DOC
+! LYP functional of the Lee, Yan, Parr, Phys. Rev B 1988, Vol 37, page 785.  
+! The expression used is the one by Miehlich, Savin, Stoll, Preuss, CPL, 1989 which gets rid of the laplacian of the density
+ END_DOC
+
+ include 'constants.include.F'
+
+! Input variables
+
+ double precision, intent(in) :: rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_2
+
+! Local variables
+
+ double precision :: a,b,c,d,c_f,omega,delta
+ double precision :: rho_13,rho_inv_13,rho_83,rho_113,rho_inv_113,denom
+ double precision :: thr,huge_num,rho_inv
+ double precision :: cst_2_113,cst_8_3,rho_2,rho_a_2,rho_b_2
+ double precision :: tmp1,tmp2,tmp3,tmp4
+ double precision :: big1,big2,big3
+
+! Output variables
+
+
+! Constants of the LYP correlation functional
+
+ a = 0.04918d0
+ b = 0.132d0
+ c = 0.2533d0
+ d = 0.349d0
+
+ thr = 1d-10
+ huge_num = 1.d0/thr
+
+ if(rho.lt.0.d0)then
+  print*,'pb !! rho.lt.0.d0'
+  stop
+ endif
+ rho_13  = rho**(1d0/3d0)
+ rho_113 = rho**(11d0/3d0)
+
+ if(dabs(rho_13) < thr) then
+  rho_inv_13 = huge_num
+ else
+  rho_inv_13 = 1.d0/rho_13
+ endif
+
+ if (dabs(rho_113) < thr) then
+  rho_inv_113 = huge_num
+ else
+  rho_inv_113 = 1d0/rho_113
+ endif
+
+ if (dabs(rho) < thr) then
+  rho_inv = huge_num
+ else
+  rho_inv = 1d0/rho
+ endif
+
+! Useful quantities to predefine
+
+ denom = 1d0/(1d0 + d*rho_inv_13)
+ omega = rho_inv_113*exp(-c*rho_inv_13)*denom
+ delta = c*rho_inv_13 + d*rho_inv_13*denom
+ c_f   = 0.3d0*(3d0*pi*pi)**(2d0/3d0)
+
+ rho_2   = rho  *rho
+ rho_a_2 = rho_a*rho_a
+ rho_b_2 = rho_b*rho_b
+
+ cst_2_113 = 2d0**(11d0/3d0)
+ cst_8_3   = 8d0/3d0
+
+ ! first term in the equation (2) of Preuss CPL, 1989
+
+ big1 = 4d0*denom*rho_a*rho_b*rho_inv
+
+ tmp1 = cst_2_113*c_f*(rho_a**cst_8_3 + rho_b**cst_8_3)
+ tmp2 = (47d0/18d0 - 7d0/18d0*delta)*grad_rho_2
+ tmp3 = - (5d0/2d0 - 1.d0/18d0*delta)*(grad_rho_a_2 + grad_rho_b_2)
+ tmp4 = - (delta - 11d0)/9d0*(rho_a*rho_inv*grad_rho_a_2 + rho_b*rho_inv*grad_rho_b_2)
+ big2 = rho_a*rho_b*(tmp1 + tmp2 + tmp3 + tmp4)
+
+ tmp1 = -2d0/3d0*rho_2*grad_rho_2
+ tmp2 = grad_rho_b_2*(2d0/3d0*rho_2 - rho_a_2)
+ tmp3 = grad_rho_a_2*(2d0/3d0*rho_2 - rho_b_2)
+ big3 = tmp1 + tmp2 + tmp3
+
+
+ ec_lyp_88 = -a*big1 -a*b*omega*big2 -a*b*omega*big3
+
+end
+
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