diff --git a/src/density.irp.f b/src/density.irp.f
index 742b768..e152ac1 100644
--- a/src/density.irp.f
+++ b/src/density.irp.f
@@ -98,4 +98,56 @@ END_PROVIDER
 
 
 
+BEGIN_PROVIDER [ double precision, density_alpha_lapl_p ]
+ implicit none
+ BEGIN_DOC
+! Laplacian of the density at the current point 
+ END_DOC
+
+ integer :: i, l
+
+ density_alpha_lapl_p = 0.
+
+ do i=1,elec_alpha_num
+  do l=1,3
+   density_alpha_lapl_p = density_alpha_lapl_p + mo_grad_p(i,l)**2
+  enddo
+  density_alpha_lapl_p = density_alpha_lapl_p + mo_value_p(i)*mo_lapl_p(i)
+ enddo
+ density_alpha_lapl_p = 2.*density_alpha_lapl_p
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, density_beta_lapl_p ]
+ implicit none
+ BEGIN_DOC
+! Laplacian of the density at the current point 
+ END_DOC
+
+ integer :: i, l
+
+ density_beta_lapl_p = 0.
+
+ do i=1,elec_beta_num
+  do l=1,3
+   density_beta_lapl_p = density_beta_lapl_p + mo_grad_p(i,l)**2
+  enddo
+  density_beta_lapl_p = density_beta_lapl_p + mo_value_p(i)*mo_lapl_p(i)
+ enddo
+ density_beta_lapl_p = 2.*density_beta_lapl_p
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, density_lapl_p ]
+ implicit none
+ BEGIN_DOC
+! Laplacian of the density at the current point 
+ END_DOC
+
+ integer :: i, l
+
+ density_lapl_p = density_alpha_lapl_p + density_beta_lapl_p
+
+END_PROVIDER