diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f
index c415ee2..0efb370 100644
--- a/src/SAMPLING/pdmc_step.irp.f
+++ b/src/SAMPLING/pdmc_step.irp.f
@@ -44,8 +44,7 @@ END_SHELL
 
   real, allocatable :: elec_coord_tmp(:,:,:)
   integer :: mod_align
-  double precision :: E_loc_save(walk_num_dmc_max)
-  double precision :: E_loc_save_tmp(walk_num_dmc_max)
+  double precision :: E_loc_save(4,walk_num_dmc_max)
   double precision :: psi_value_save(walk_num)
   double precision :: psi_value_save_tmp(walk_num)
   double precision :: pdmc_weight(walk_num)
@@ -54,7 +53,6 @@ END_SHELL
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save_tmp
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: E_loc_save
-  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: E_loc_save_tmp
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save_tmp
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: pdmc_weight
@@ -131,7 +129,7 @@ END_SHELL
           psi_grad_psi_inv_z(i) = psi_grad_psi_inv_save(i,3,i_walk)
         enddo
         psi_value = psi_value_save(i_walk)
-        E_loc = E_loc_save(i_walk)
+        E_loc = E_loc_save(1,i_walk)
       enddo
       SOFT_TOUCH elec_coord psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z psi_value E_loc
     else
@@ -142,7 +140,7 @@ END_SHELL
       enddo
       TOUCH elec_coord
       psi_value_save(i_walk) = psi_value
-      E_loc_save(i_walk) = E_loc
+      E_loc_save(:,i_walk) = E_loc
     endif
 
    double precision               :: p,q
@@ -151,13 +149,19 @@ END_SHELL
    call brownian_step(p,q,accepted,delta_x)
 
 !   if ( psi_value * psi_value_save(i_walk) >= 0.d0 ) then
-     delta = ((E_loc+E_loc_save(i_walk))*0.5d0 - E_ref) * p
-     if ( delta > thr ) then
-       pdmc_weight(i_walk) = dexp(-dtime_step*thr)
-     else if ( delta < -thr ) then
-       pdmc_weight(i_walk) = dexp(dtime_step*thr)
-     else
+
+!2    delta = (E_loc+E_loc_save(1,i_walk))*0.5d0
+!3    delta = (5.d0 * E_loc + 8.d0 * E_loc_save(1,i_walk) - E_loc_save(2,i_walk))/12.d0
+      delta = (9.d0*E_loc+19.d0*E_loc_save(1,i_walk)-5.d0*E_loc_save(2,i_walk)+E_loc_save(3,i_walk))/24.d0
+!      delta = -((-251.d0*E_loc)-646.d0*E_loc_save(1,i_walk)+264.d0*E_loc_save(2,i_walk)-&
+!         106.d0*E_loc_save(3,i_walk)+19.d0*E_loc_save(4,i_walk))/720.d0
+
+     delta = (delta - E_ref)*p
+       
+     if (delta >= 0.d0) then
        pdmc_weight(i_walk) = dexp(-dtime_step*delta)
+     else
+       pdmc_weight(i_walk) = 2.d0-dexp(dtime_step*delta)
      endif
      elec_coord(elec_num+1,1) += p*time_step
      elec_coord(elec_num+1,2)  = E_loc
@@ -174,7 +178,10 @@ END_SHELL
      enddo
 
      psi_value_save(i_walk) = psi_value
-     E_loc_save(i_walk) = E_loc
+     E_loc_save(4,i_walk) = E_loc_save(3,i_walk)
+     E_loc_save(3,i_walk) = E_loc_save(2,i_walk)
+     E_loc_save(2,i_walk) = E_loc_save(1,i_walk)
+     E_loc_save(1,i_walk) = E_loc
 
      if (dabs(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag)) > 1.d-15) then
        dmc_zv_weight = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag))