diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f
index 2764609..7044080 100644
--- a/src/SAMPLING/dmc_step.irp.f
+++ b/src/SAMPLING/dmc_step.irp.f
@@ -248,14 +248,15 @@ END_SHELL
   if (cpu1 < cpu0) then
     cpu1 = cpu1+cpu0
   endif
-  loop = dble(cpu1-cpu0) < dble(block_time)*dble(count_rate)
+  loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time
   if (cpu1-cpu2 > count_rate) then
     integer                        :: do_run
     call get_running(do_run)
-    loop = do_run == t_Running
+    loop = loop.and.(do_run == t_Running)
     cpu2 = cpu1
   endif
 
+
   SOFT_TOUCH elec_coord_full_dmc psi_value psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z elec_coord 
 
  enddo
diff --git a/src/SAMPLING/fkmc_step.irp.f b/src/SAMPLING/fkmc_step.irp.f
index 878af69..3f57d20 100644
--- a/src/SAMPLING/fkmc_step.irp.f
+++ b/src/SAMPLING/fkmc_step.irp.f
@@ -319,11 +319,11 @@ END_SHELL
   if (cpu1 < cpu0) then
     cpu1 = cpu1+cpu0
   endif
-  loop = dble(cpu1-cpu0) < dble(block_time)*dble(count_rate)
+  loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time
   if (cpu1-cpu2 > count_rate) then
     integer                        :: do_run
     call get_running(do_run)
-    loop = do_run == t_Running
+    loop = loop.and.(do_run == t_Running)
     cpu2 = cpu1
   endif
 
diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f
new file mode 100644
index 0000000..e4c3e70
--- /dev/null
+++ b/src/SAMPLING/pdmc_step.irp.f
@@ -0,0 +1,380 @@
+! Providers of *_pdmc_block_walk
+!==============================
+BEGIN_SHELL [ /usr/bin/python ]
+from properties import *
+
+t = """
+ BEGIN_PROVIDER [ $T, $X_pdmc_block_walk $D1 ]
+&BEGIN_PROVIDER [ $T, $X_pdmc_block_walk_kahan $D2 ]
+&BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk $D1 ]
+&BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk_kahan $D2 ]
+ implicit none
+ BEGIN_DOC  
+! pdMC averages of $X. Computed in E_loc_pdmc_block_walk
+ END_DOC
+ $X_pdmc_block_walk = 0.d0
+ $X_pdmc_block_walk_kahan = 0.d0
+ $X_2_pdmc_block_walk = 0.d0
+ $X_2_pdmc_block_walk_kahan = 0.d0
+END_PROVIDER
+"""
+for p in properties:
+  if p[1] != 'e_loc':
+    if p[2] == "":
+      D1 = ""
+      D2 = ", (3)"
+    else:
+      D1 = ", ("+p[2][1:-1]+")"
+      D2 = ", ("+p[2][1:-1]+",3)"
+    print t.replace("$X",p[1]).replace("$T",p[0]).replace("$D1",D1).replace("$D2",D2)
+
+END_SHELL
+
+
+
+ BEGIN_PROVIDER [ double precision, E_loc_pdmc_block_walk       ]
+&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk     ]
+&BEGIN_PROVIDER [ double precision, E_loc_pdmc_block_walk_kahan, (3) ]
+&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk_kahan, (3) ]
+ implicit none
+ include '../types.F'
+ BEGIN_DOC
+! Properties averaged over the block using the PDMC method
+ END_DOC
+
+  real, allocatable :: elec_coord_tmp(:,:,:)
+  integer :: mod_align
+  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)
+  double precision, allocatable :: psi_grad_psi_inv_save(:,:,:)
+  double precision, allocatable :: psi_grad_psi_inv_save_tmp(:,:,:)
+  !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 :: psi_value_save
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save_tmp
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: pdmc_weight
+  allocate ( psi_grad_psi_inv_save(elec_num_8,3,walk_num) ,          &
+       psi_grad_psi_inv_save_tmp(elec_num_8,3,walk_num) ,            &
+      elec_coord_tmp(mod_align(elec_num+1),3,walk_num) )
+  psi_value_save = 0.d0
+  psi_value_save_tmp = 0.d0
+  pdmc_weight = 1.d0
+
+! Initialization
+ if (vmc_algo /= t_Brownian) then
+   call abrt(irp_here,'PDMC should run with Brownian algorithm')
+ endif
+
+ integer :: k, i_walk, i_step
+
+BEGIN_SHELL [ /usr/bin/python ]
+from properties import *
+t = """
+ if (calc_$X) then
+   !DIR$ VECTOR ALIGNED
+   $X_pdmc_block_walk = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_pdmc_block_walk_kahan = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_2_pdmc_block_walk = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_2_pdmc_block_walk_kahan = 0.d0
+ endif
+"""
+for p in properties:
+ print  t.replace("$X",p[1])
+END_SHELL
+
+ logical                        :: loop
+ integer*8                      :: cpu0, cpu1, cpu2, count_rate, count_max
+
+ loop = .True.
+ call system_clock(cpu0, count_rate, count_max)
+ cpu2 = cpu0
+
+ block_weight = 0.d0
+
+ real, external                 :: accep_rate
+ double precision               :: delta, thr
+
+ thr = 2.d0/time_step_sq
+
+ logical :: first_loop
+ first_loop = .True.
+ if (walk_num > 1) then
+   call abrt(irp_here,'walk_num > 1')
+ endif
+
+  integer :: info
+!  double precision :: H(0:pdmc_n_diag/2,0:pdmc_n_diag/2), S(0:pdmc_n_diag/2,0:pdmc_n_diag/2), w(0:pdmc_n_diag/2), work(3*pdmc_n_diag+1) 
+!  H = 0.d0
+!  S = 0.d0
+
+ do while (loop)
+
+  i_walk = 1
+    
+    if (.not.first_loop) then
+      integer                        :: i,j,l
+      do l=1,3
+        do i=1,elec_num+1
+          elec_coord(i,l) = elec_coord_full(i,l,i_walk)
+        enddo
+        do i=1,elec_num
+          psi_grad_psi_inv_x(i) = psi_grad_psi_inv_save(i,1,i_walk)
+          psi_grad_psi_inv_y(i) = psi_grad_psi_inv_save(i,2,i_walk)
+          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(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
+      do l=1,3
+        do i=1,elec_num+1
+          elec_coord(i,l) = elec_coord_full(i,l,i_walk)
+        enddo
+      enddo
+      TOUCH elec_coord
+      psi_value_save(i_walk) = psi_value
+      E_loc_save(:,i_walk) = E_loc
+    endif
+
+   double precision               :: p,q
+   real                           :: delta_x
+   logical                        :: accepted
+   call brownian_step(p,q,accepted,delta_x)
+
+!   if ( psi_value * psi_value_save(i_walk) >= 0.d0 ) then
+
+!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
+     elec_coord(elec_num+1,3)  = pdmc_weight(i_walk) * pdmc_pop_weight_mult(pdmc_n_diag)
+     do l=1,3
+        do i=1,elec_num+1
+          elec_coord_full(i,l,i_walk) = elec_coord(i,l)
+        enddo
+     enddo
+     do i=1,elec_num
+       psi_grad_psi_inv_save(i,1,i_walk) = psi_grad_psi_inv_x(i)
+       psi_grad_psi_inv_save(i,2,i_walk) = psi_grad_psi_inv_y(i)
+       psi_grad_psi_inv_save(i,3,i_walk) = psi_grad_psi_inv_z(i)
+     enddo
+
+     psi_value_save(i_walk) = psi_value
+     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))
+       dmc_zv_weight_half = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag/2))
+     else
+       dmc_zv_weight = 0.d0
+       dmc_zv_weight_half = 0.d0
+     endif
+     TOUCH dmc_zv_weight dmc_zv_weight_half
+
+!    do i=1,pdmc_n_diag+1
+!      E_loc_zv(i) =   E_loc * pdmc_pop_weight_mult(i-1) * pdmc_weight(i_walk) * dmc_zv_weight + (E_trial-E_loc) * dmc_zv_weight
+!      E_loc_zv(i+pdmc_n_diag+1) = pdmc_pop_weight_mult(i-1) * pdmc_weight(i_walk) * dmc_zv_weight
+!    enddo
+
+BEGIN_SHELL [ /usr/bin/python ]
+from properties import *
+t = """
+     if (calc_$X) then
+   ! Kahan's summation algorithm to compute these sums reducing the rounding error:
+   !  $X_pdmc_block_walk    += $X * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk)
+   !  $X_2_pdmc_block_walk  += $X_2 * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk)
+   ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
+   
+      $X_pdmc_block_walk_kahan($D2 3) = $X * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) - $X_pdmc_block_walk_kahan($D2 1)
+      $X_pdmc_block_walk_kahan($D2 2) = $X_pdmc_block_walk $D1  + $X_pdmc_block_walk_kahan($D2 3)
+      $X_pdmc_block_walk_kahan($D2 1) = ($X_pdmc_block_walk_kahan($D2 2) - $X_pdmc_block_walk $D1 ) &
+          - $X_pdmc_block_walk_kahan($D2 3)
+      $X_pdmc_block_walk $D1  =  $X_pdmc_block_walk_kahan($D2 2) 
+   
+   
+      $X_2_pdmc_block_walk_kahan($D2 3) = $X_2 * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) - $X_2_pdmc_block_walk_kahan($D2 1)
+      $X_2_pdmc_block_walk_kahan($D2 2) = $X_2_pdmc_block_walk $D1 + $X_2_pdmc_block_walk_kahan($D2 3)
+      $X_2_pdmc_block_walk_kahan($D2 1) = ($X_2_pdmc_block_walk_kahan($D2 2) - $X_2_pdmc_block_walk $D1 ) &
+          - $X_2_pdmc_block_walk_kahan($D2 3)
+      $X_2_pdmc_block_walk $D1 =  $X_2_pdmc_block_walk_kahan($D2 2) 
+     endif
+"""
+for p in properties:
+  if p[2] == "":
+   D1 = ""
+   D2 = ""
+  else:
+   D1 = "("+":"*(p[2].count(',')+1)+")"
+   D2 = ":"*(p[2].count(',')+1)+","
+  print t.replace("$X",p[1]).replace("$D1",D1).replace("$D2",D2)
+
+END_SHELL
+
+    block_weight += pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk)
+
+    pdmc_pop_weight_mult(0) = 1.d0/pdmc_weight(i_walk)
+!    do k=0,pdmc_n_diag/2
+!      do l=0,pdmc_n_diag/2
+!        H(k,l) += E_loc*pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk)
+!        S(k,l) += pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk)
+!      enddo
+!    enddo
+!    H = H + (E_trial - E_loc)
+
+!  else
+!     pdmc_weight(i_walk) = 1.d0
+!     pdmc_pop_weight(:,:) = 1.d0
+!     pdmc_pop_weight_mult(:) = 1.d0
+!  endif
+   
+
+  do k=1,pdmc_n_diag
+    ! Move to the next projection step
+    if (pdmc_projection(pdmc_n_diag) > 0) then
+      pdmc_projection_step(k) = mod(pdmc_projection_step(k),pdmc_projection(k))+1
+    else
+      pdmc_projection_step(k) = 1
+    endif
+
+    ! Eventually, recompute the weight of the population
+    if (pdmc_projection_step(k) == k) then
+      pdmc_pop_weight_mult(k) = 1.d0
+      do l=1,pdmc_projection(k)
+        pdmc_pop_weight_mult(k) *= pdmc_pop_weight(l,k)
+      enddo
+    endif
+    ! Remove contribution of the old value of the weight at the new
+    ! projection step
+    
+    pdmc_pop_weight_mult(k) *= 1.d0/pdmc_pop_weight(pdmc_projection_step(k),k)
+    pdmc_pop_weight(pdmc_projection_step(k),k) = pdmc_weight(i_walk)/dble(walk_num)
+
+    ! Update the running population weight
+    pdmc_pop_weight_mult(k) *= pdmc_pop_weight(pdmc_projection_step(k),k)
+
+  enddo
+
+
+  call system_clock(cpu1, count_rate, count_max)
+  if (cpu1 < cpu0) then
+    cpu1 = cpu1+cpu0
+  endif
+  loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time
+  if (cpu1-cpu2 > count_rate) then
+    integer                        :: do_run
+    call get_running(do_run)
+    loop = loop.and.(do_run == t_Running)
+    cpu2 = cpu1
+  endif
+
+  SOFT_TOUCH elec_coord_full pdmc_pop_weight_mult 
+
+  first_loop = .False.
+
+ enddo
+
+ 
+ double precision :: factor
+ factor = 1.d0/block_weight
+ SOFT_TOUCH block_weight
+
+BEGIN_SHELL [ /usr/bin/python ]
+from properties import *
+t = """
+ if (calc_$X) then
+   $X_pdmc_block_walk   *= factor
+   $X_2_pdmc_block_walk *= factor
+ endif
+"""
+for p in properties:
+ print  t.replace("$X",p[1])
+END_SHELL
+
+!  H(0,0) = H(3,3)
+!  H(1,0) = H(4,3)
+!  H(0,1) = H(3,4)
+!  H(1,1) = H(4,4)
+!  S(0,0) = S(3,3)
+!  S(1,0) = S(4,3)
+!  S(0,1) = S(3,4)
+!  S(1,1) = S(4,4)
+!
+!  print *,  H(0,0)/S(0,0)
+!  print *,  H(1,1)/S(1,1)
+!  print *,  ''
+!
+!  call dsygv(1, 'N', 'U', pdmc_n_diag/2+1, H, pdmc_n_diag/2+1, S, pdmc_n_diag/2+1, w, work, 3*(pdmc_n_diag+1), info)
+!  call dsygv(1, 'N', 'U', 2, H, pdmc_n_diag/2+1, S, pdmc_n_diag/2+1, w, work, 3*(pdmc_n_diag+1), info)
+!  E_loc_zv_diag_pdmc_block_walk = w(0)
+!  print *,  w
+
+ deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp )
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer, pdmc_projection, (pdmc_n_diag) ]
+&BEGIN_PROVIDER [ integer, pdmc_projection_step, (pdmc_n_diag) ]
+ implicit none
+ BEGIN_DOC  
+! Number of projection steps for PDMC
+ END_DOC
+ real :: pdmc_projection_time
+ pdmc_projection_time = 1.
+ call get_simulation_srmc_projection_time(pdmc_projection_time)
+ pdmc_projection(pdmc_n_diag) = int( pdmc_projection_time/time_step)
+ integer :: k
+ do k=1,pdmc_n_diag-1
+   pdmc_projection(k) = k*pdmc_projection(pdmc_n_diag)/pdmc_n_diag
+ enddo
+ pdmc_projection_step(:) = 0
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection(pdmc_n_diag)+1,pdmc_n_diag) ]
+ implicit none
+ BEGIN_DOC  
+! Population weight of PDMC
+ END_DOC
+ pdmc_pop_weight(:,:) = 1.d0
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, pdmc_pop_weight_mult, (0:pdmc_n_diag) ]
+ implicit none
+ BEGIN_DOC  
+! Population weight of PDMC
+ END_DOC
+ pdmc_pop_weight_mult(:) = 1.d0
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, pdmc_n_diag ]
+ implicit none
+ BEGIN_DOC
+! Size of the matrix to diagonalize
+ END_DOC
+ pdmc_n_diag = 8
+END_PROVIDER
+
+
+
diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f
index 836ddae..1300fe2 100644
--- a/src/SAMPLING/srmc_step.irp.f
+++ b/src/SAMPLING/srmc_step.irp.f
@@ -319,11 +319,11 @@ END_SHELL
   if (cpu1 < cpu0) then
     cpu1 = cpu1+cpu0
   endif
-  loop = dble(cpu1-cpu0) < dble(block_time)*dble(count_rate)
+  loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time
   if (cpu1-cpu2 > count_rate) then
     integer                        :: do_run
     call get_running(do_run)
-    loop = do_run == t_Running
+    loop = loop.and.(do_run == t_Running)
     cpu2 = cpu1
   endif
 
diff --git a/src/SAMPLING/vmc_step.irp.f b/src/SAMPLING/vmc_step.irp.f
index 7aedf42..e09ee80 100644
--- a/src/SAMPLING/vmc_step.irp.f
+++ b/src/SAMPLING/vmc_step.irp.f
@@ -132,11 +132,11 @@ END_SHELL
       if (cpu1 < cpu0) then
         cpu1 = cpu1+cpu0
       endif
-      loop = dble(cpu1-cpu0)*dble(walk_num) < dble(block_time)*dble(count_rate)
+      loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time
       if (cpu1-cpu2 > count_rate) then
         integer                        :: do_run
         call get_running(do_run)
-        loop = do_run == t_Running
+        loop = loop.and.(do_run == t_Running)
         cpu2 = cpu1
       endif