diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f
index 12ebb73..56dc6ab 100644
--- a/src/PROPERTIES/properties_energy.irp.f
+++ b/src/PROPERTIES/properties_energy.irp.f
@@ -13,7 +13,7 @@ BEGIN_PROVIDER [ double precision, single_det_E_kin ]
   do i=1,elec_num
     single_det_E_kin -= 0.5d0*single_det_lapl(i)/single_det_value
   enddo
-  
+
 END_PROVIDER
 
 
@@ -22,7 +22,7 @@ BEGIN_PROVIDER  [ double precision, single_det_E_loc ]
   BEGIN_DOC
   ! Local energy : single_det_E_kin + E_pot + E_nucl
   END_DOC
-  
+
   single_det_E_loc = single_det_E_kin + E_pot + E_nucl
 END_PROVIDER
 
@@ -32,7 +32,7 @@ BEGIN_PROVIDER [ double precision, E_pot_grad, (elec_num,3) ]
   BEGIN_DOC
   ! Gradient of the Electronic Potential energy
   END_DOC
-  
+
   integer                        :: i,j
   double precision               :: dinv
   do i=1,elec_num
@@ -64,7 +64,7 @@ BEGIN_PROVIDER [ double precision, E_pot_grad, (elec_num,3) ]
       E_pot_grad(i,3) += nucl_elec_dist_vec(3,j,i)*dinv
     enddo
   enddo
-  
+
 END_PROVIDER
 
 
@@ -73,7 +73,7 @@ BEGIN_PROVIDER [ double precision, E_pot_elec, (elec_num) ]
   BEGIN_DOC
   ! Electronic Potential energy
   END_DOC
-  
+
   integer                        :: i, j
   if (do_pseudo) then
     do i=1,elec_num
@@ -89,7 +89,7 @@ BEGIN_PROVIDER [ double precision, E_pot_elec, (elec_num) ]
     !DIR$ VECTOR ALIGNED
     !DIR$ LOOP COUNT(50)
     do j=1,elec_num
-      E_pot_elec(i) = E_pot_elec(i) + 0.5d0*elec_dist_inv(j,i) 
+      E_pot_elec(i) = E_pot_elec(i) + 0.5d0*elec_dist_inv(j,i)
     enddo
     !DIR$ VECTOR ALIGNED
     !DIR$ LOOP COUNT(50)
@@ -97,7 +97,7 @@ BEGIN_PROVIDER [ double precision, E_pot_elec, (elec_num) ]
       E_pot_elec(i) = E_pot_elec(i) - nucl_charge(j)*nucl_elec_dist_inv(j,i)
     enddo
   enddo
-  
+
 END_PROVIDER
 
 
@@ -106,7 +106,7 @@ BEGIN_PROVIDER [ double precision, E_pot_elec_one, (elec_num) ]
   BEGIN_DOC
   ! Electronic Potential energy
   END_DOC
-  
+
   integer                        :: i, j
   do i=1,elec_num
     E_pot_elec_one(i) = 0.d0
@@ -116,7 +116,7 @@ BEGIN_PROVIDER [ double precision, E_pot_elec_one, (elec_num) ]
       E_pot_elec_one(i) -= nucl_charge(j)*nucl_elec_dist_inv(j,i)
     enddo
   enddo
-  
+
 END_PROVIDER
 
 
@@ -125,7 +125,7 @@ BEGIN_PROVIDER [ double precision, E_pot_elec_two, (elec_num) ]
   BEGIN_DOC
   ! Electronic Potential energy
   END_DOC
-  
+
   integer                        :: i, j
   do i=1,elec_num
     E_pot_elec_two(i) = 0.d0
@@ -138,22 +138,22 @@ BEGIN_PROVIDER [ double precision, E_pot_elec_two, (elec_num) ]
       E_pot_elec_two(i) += 0.5d0*elec_dist_inv(j,i)
     enddo
   enddo
-  
+
 END_PROVIDER
 
 
 BEGIN_PROVIDER [ double precision, E_kin_elec, (elec_num) ]
   implicit none
-  
+
   BEGIN_DOC
   ! Electronic Kinetic energy : -1/2 (Lapl.Psi)/Psi
   END_DOC
-  
+
   integer                        :: i
   do i=1,elec_num
     E_kin_elec(i) = -0.5d0*psi_lapl_psi_inv(i)
   enddo
-  
+
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, dmc_zv_weight ]
@@ -182,7 +182,7 @@ BEGIN_PROVIDER [ double precision, E_nucl ]
   BEGIN_DOC
   ! Nuclear potential energy
   END_DOC
-  
+
   E_nucl = 0.d0
   integer                        :: i, j
   do i=1,nucl_num
@@ -190,7 +190,7 @@ BEGIN_PROVIDER [ double precision, E_nucl ]
       E_nucl += nucl_charge(i)*nucl_charge(j)/nucl_dist(j,i)
     enddo
   enddo
-  
+
   E_nucl_min = min(E_nucl,E_nucl_min)
   E_nucl_max = max(E_nucl,E_nucl_max)
   SOFT_TOUCH E_nucl_min E_nucl_max
@@ -202,13 +202,13 @@ BEGIN_PROVIDER [ double precision, E_pot ]
   BEGIN_DOC
   ! Electronic Potential energy
   END_DOC
-  
+
   E_pot = 0.d0
   integer                        :: i, j
   do i=1,elec_num
     E_pot += E_pot_elec(i)
   enddo
-  
+
   E_pot_min = min(E_pot,E_pot_min)
   E_pot_max = max(E_pot,E_pot_max)
   SOFT_TOUCH E_pot_min E_pot_max
@@ -220,16 +220,16 @@ BEGIN_PROVIDER [ double precision, E_kin ]
   BEGIN_DOC
   ! Electronic Kinetic energy : -1/2 (Lapl.Psi)/Psi
   END_DOC
-  
+
   E_kin = 0.d0
-  
+
   integer                        :: i
   !DIR$ VECTOR ALIGNED
   !DIR$ LOOP COUNT(200)
   do i=1,elec_num
     E_kin -= 0.5d0*psi_lapl_psi_inv(i)
   enddo
-  
+
   E_kin_min = min(E_kin,E_kin_min)
   E_kin_max = max(E_kin,E_kin_max)
   SOFT_TOUCH E_kin_min E_kin_max
@@ -242,7 +242,7 @@ BEGIN_PROVIDER  [ double precision, E_loc ]
   BEGIN_DOC
   ! Local energy : E_kin + E_pot + E_nucl
   END_DOC
-  
+
   integer                        :: i
   E_loc = E_nucl
   !DIR$ VECTOR ALIGNED
@@ -250,14 +250,14 @@ BEGIN_PROVIDER  [ double precision, E_loc ]
   do i=1,elec_num
     E_loc += E_kin_elec(i) + E_pot_elec(i)
   enddo
-  
-!  ! Avoid divergence of E_loc
-!  if (qmc_method == t_DMC) then
-!    double precision :: delta_e
+
+  ! Avoid divergence of E_loc and population explosion
+  if (do_pseudo) then
+    double precision :: delta_e
 !    delta_e = E_loc-E_ref
-!    E_loc = E_ref + erf(1.d0/(time_step*delta_e*time_step*delta_e)) * delta_e
-!  endif
-  
+!    E_loc = E_ref + erf(delta_e*time_step_sq)/time_step_sq
+     E_loc = max(2.d0*E_ref, E_loc)
+  endif
   E_loc_min = min(E_loc,E_loc_min)
   E_loc_max = max(E_loc,E_loc_max)
   SOFT_TOUCH E_loc_min E_loc_max
@@ -270,7 +270,7 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ]
  BEGIN_DOC
  ! Zero-variance parameter on E_loc
  END_DOC
-   E_loc_zv = E_loc 
+   E_loc_zv = E_loc
    E_loc_zv += (E_trial-E_loc) * dmc_zv_weight
 !   E_loc_zv += - time_step*(E_trial**2 + 1.44341217940434 - E_loc**2)*dmc_zv_weight
 !  E_loc_zv(3)  = dmc_zv_weight_half