diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f
new file mode 100644
index 0000000..4cd823f
--- /dev/null
+++ b/src/SAMPLING/srmc_step.irp.f
@@ -0,0 +1,306 @@
+! Providers of *_srmc_block_walk
+!==============================
+BEGIN_SHELL [ /usr/bin/python ]
+from properties import *
+
+t = """
+ BEGIN_PROVIDER [ $T, $X_srmc_block_walk $D1 ]
+&BEGIN_PROVIDER [ $T, $X_srmc_block_walk_kahan $D2 ]
+&BEGIN_PROVIDER [ $T, $X_2_srmc_block_walk $D1 ]
+&BEGIN_PROVIDER [ $T, $X_2_srmc_block_walk_kahan $D2 ]
+ implicit none
+ BEGIN_DOC  
+! SRMC averages of $X. Computed in E_loc_srmc_block_walk
+ END_DOC
+ $X_srmc_block_walk = 0.d0
+ $X_srmc_block_walk_kahan = 0.d0
+ $X_2_srmc_block_walk = 0.d0
+ $X_2_srmc_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_srmc_block_walk       ]
+&BEGIN_PROVIDER [ double precision, E_loc_2_srmc_block_walk     ]
+&BEGIN_PROVIDER [ double precision, E_loc_srmc_block_walk_kahan, (3) ]
+&BEGIN_PROVIDER [ double precision, E_loc_2_srmc_block_walk_kahan, (3) ]
+ implicit none
+ include '../types.F'
+ BEGIN_DOC
+! Properties averaged over the block using the SRMC method
+ END_DOC
+
+  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 :: psi_value_save(walk_num)
+  double precision :: psi_value_save_tmp(walk_num)
+  double precision :: srmc_weight(walk_num)
+  !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 :: srmc_weight
+  allocate ( elec_coord_tmp(mod_align(elec_num+1),3,walk_num) )
+  psi_value_save = 0.d0
+  psi_value_save_tmp = 0.d0
+  srmc_weight = 1.d0
+
+! Initialization
+ if (vmc_algo /= t_Brownian) then
+   call abrt(irp_here,'SRMC 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_srmc_block_walk = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_srmc_block_walk_kahan = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_2_srmc_block_walk = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_2_srmc_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
+
+
+ do while (loop)
+
+  ! Every walker makes a step
+  do i_walk=1,walk_num
+    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
+    enddo
+    TOUCH elec_coord
+
+
+   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
+     delta = ((E_loc+E_loc_save(i_walk))*0.5d0 - E_ref) * p
+     if ( delta > thr ) then
+       srmc_weight(i_walk) = dexp(-dtime_step*thr)
+     else if ( delta < -thr ) then
+       srmc_weight(i_walk) = dexp(dtime_step*thr)
+     else
+       srmc_weight(i_walk) = dexp(-dtime_step*delta)
+     endif
+   else
+     srmc_weight(i_walk) = 0.d0
+   endif
+   
+   elec_coord(elec_num+1,1) += p*time_step
+   elec_coord(elec_num+1,2)  = E_loc
+   elec_coord(elec_num+1,3)  = srmc_weight(i_walk)
+   do l=1,3
+      do i=1,elec_num+1
+        elec_coord_full(i,l,i_walk) = elec_coord(i,l)
+      enddo
+   enddo
+
+   psi_value_save(i_walk) = psi_value
+   E_loc_save(i_walk) = E_loc
+
+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_srmc_block_walk    += $X * pop_weight_mult * srmc_weight(i_walk)
+   !  $X_2_srmc_block_walk  += $X_2 * pop_weight_mult * srmc_weight(i_walk)
+   ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
+   
+      $X_srmc_block_walk_kahan($D2 3) = $X * pop_weight_mult * srmc_weight(i_walk) - $X_srmc_block_walk_kahan($D2 1)
+      $X_srmc_block_walk_kahan($D2 2) = $X_srmc_block_walk $D1  + $X_srmc_block_walk_kahan($D2 3)
+      $X_srmc_block_walk_kahan($D2 1) = ($X_srmc_block_walk_kahan($D2 2) - $X_srmc_block_walk $D1 ) &
+          - $X_srmc_block_walk_kahan($D2 3)
+      $X_srmc_block_walk $D1  =  $X_srmc_block_walk_kahan($D2 2) 
+   
+   
+      $X_2_srmc_block_walk_kahan($D2 3) = $X_2 * pop_weight_mult * srmc_weight(i_walk) - $X_2_srmc_block_walk_kahan($D2 1)
+      $X_2_srmc_block_walk_kahan($D2 2) = $X_2_srmc_block_walk $D1 + $X_2_srmc_block_walk_kahan($D2 3)
+      $X_2_srmc_block_walk_kahan($D2 1) = ($X_2_srmc_block_walk_kahan($D2 2) - $X_2_srmc_block_walk $D1 ) &
+          - $X_2_srmc_block_walk_kahan($D2 3)
+      $X_2_srmc_block_walk $D1 =  $X_2_srmc_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 += pop_weight_mult * srmc_weight(i_walk)
+
+  enddo
+
+  ! Move to the next projection step
+  if (srmc_projection > 0) then
+    srmc_projection_step = mod(srmc_projection_step,srmc_projection)+1
+  else
+    srmc_projection_step = 1
+  endif
+
+  ! Eventually, recompute the weight of the population
+  if (srmc_projection_step == 1) then
+    pop_weight_mult = 1.d0
+    do k=1,srmc_projection
+      pop_weight_mult *= pop_weight(k)
+    enddo
+  endif
+
+  ! Remove contribution of the old value of the weight at the new
+  ! projection step
+  pop_weight_mult *= 1.d0/pop_weight(srmc_projection_step)
+
+  ! Compute the new weight of the population
+  double precision :: sum_weight
+  sum_weight = 0.d0
+  do k=1,walk_num
+    sum_weight += srmc_weight(k)
+  enddo
+  pop_weight(srmc_projection_step) = sum_weight/dble(walk_num)
+
+  ! Update the running population weight
+  pop_weight_mult *= pop_weight(srmc_projection_step)
+
+! Reconfiguration
+  integer :: ipos(walk_num)
+
+  call reconfigure(ipos,srmc_weight)
+   
+  do k=1,walk_num
+    do l=1,3
+     do i=1,elec_num+1
+      elec_coord_tmp(i,l,k) = elec_coord_full(i,l,k)
+     enddo
+    enddo
+    psi_value_save_tmp(k) = psi_value_save(k)
+    E_loc_save_tmp(k) = E_loc_save(k)
+  enddo
+  
+  integer :: ipm
+  do k=1,walk_num
+   ipm = ipos(k)
+   do l=1,3
+    do i=1,elec_num+1
+     elec_coord_full(i,l,k) = elec_coord_tmp(i,l,ipm)
+    enddo
+   enddo
+   psi_value_save(k) = psi_value_save_tmp(ipm)
+   E_loc_save(k) = E_loc_save_tmp(ipm)
+  enddo
+
+  call system_clock(cpu1, count_rate, count_max)
+  if (cpu1 < cpu0) then
+    cpu1 = cpu1+cpu0
+  endif
+  loop = dble(cpu1-cpu0) < dble(block_time)*dble(count_rate)
+  if (cpu1-cpu2 > count_rate) then
+    integer                        :: do_run
+    call get_running(do_run)
+    loop = do_run == t_Running
+    cpu2 = cpu1
+  endif
+
+  SOFT_TOUCH elec_coord_full psi_value psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z elec_coord pop_weight_mult
+
+ 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_srmc_block_walk   *= factor
+   $X_2_srmc_block_walk *= factor
+ endif
+"""
+for p in properties:
+ print  t.replace("$X",p[1])
+END_SHELL
+
+ deallocate ( elec_coord_tmp )
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, pop_weight_mult ]
+ implicit none
+ BEGIN_DOC  
+! Population weight of SRMC
+ END_DOC
+ pop_weight_mult = pop_weight(srmc_projection)
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, srmc_projection ]
+&BEGIN_PROVIDER [ integer, srmc_projection_step ]
+ implicit none
+ BEGIN_DOC  
+! Number of projection steps for SRMC
+ END_DOC
+ real :: srmc_projection_time
+ srmc_projection_time = 1.
+ call get_simulation_srmc_projection_time(srmc_projection_time)
+ srmc_projection = int( srmc_projection_time/time_step)
+ srmc_projection_step = 0
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, pop_weight, (0:srmc_projection+1) ]
+ implicit none
+ BEGIN_DOC  
+! Population weight of SRMC
+ END_DOC
+ pop_weight = 1.d0
+END_PROVIDER
+
+