Better time step error in SRMC

src/SAMPLING/srmc_step.irp.f

@@ -44,8 +44,8 @@ END_SHELL
 
   real, allocatable :: elec_coord_tmp(:,:,:)
   integer :: mod_align
-  double precision :: E_loc_save(walk_num_dmc_max)
+  double precision :: E_loc_save(4,walk_num_dmc_max)
-  double precision :: E_loc_save_tmp(walk_num_dmc_max)
+  double precision :: E_loc_save_tmp(4,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)
@@ -124,7 +124,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
@@ -135,7 +135,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
@@ -144,19 +144,54 @@ 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
+!    delta = (E_loc+E_loc_save(1,i_walk))*0.5d0
-     if ( delta > thr ) then
+!    delta = (5.d0 * E_loc + 8.d0 * E_loc_save(1,i_walk) - E_loc_save(2,i_walk))/12.d0
-       srmc_weight(i_walk) = dexp(-dtime_step*thr)
+
-     else if ( delta < -thr ) then
+    delta = (9.d0*E_loc+19.d0*E_loc_save(1,i_walk)- &
-       srmc_weight(i_walk) = dexp(dtime_step*thr)
+             5.d0*E_loc_save(2,i_walk)+E_loc_save(3,i_walk))/24.d0
-     else
+
+!      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
        srmc_weight(i_walk) = dexp(-dtime_step*delta)
+     else
+       srmc_weight(i_walk) = 2.d0-dexp(dtime_step*delta)
      endif
+
+
+!    if (accepted) then
+!      ! Compute correction to past weights
+!      double precision               :: delta_old, delta_new
+!      delta_old = (9.d0*E_loc_save(1,i_walk)+19.d0*E_loc_save(2,i_walk)-&
+!          5.d0*E_loc_save(3,i_walk)+E_loc_save(4,i_walk))/24.d0 - E_ref
+!      
+!      
+!      if (delta_old >= 0.d0) then
+!        srmc_weight(i_walk) = srmc_weight(i_walk) * dexp(dtime_step*delta_old)
+!      else
+!        srmc_weight(i_walk) = srmc_weight(i_walk) * (2.d0-dexp(-dtime_step*delta_old))
+!      endif
+!      
+!      delta_new = (-(E_loc_save_tmp(3,i_walk)-13.d0*E_loc_save_tmp(2,i_walk)&
+!          -13.d0*E_loc_save_tmp(1,i_walk)+E_loc))/24.d0  - E_ref
+!      
+!      if (delta_new >= 0.d0) then
+!        srmc_weight(i_walk) = srmc_weight(i_walk) * dexp(-dtime_step*delta_new)
+!      else
+!        srmc_weight(i_walk) = srmc_weight(i_walk) * (2.d0-dexp(dtime_step*delta_new) )
+!      endif
+!      
+!    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) * pop_weight_mult
      do l=1,3
         do i=1,elec_num+1
+
           elec_coord_full(i,l,i_walk) = elec_coord(i,l)
         enddo
      enddo
@@ -167,25 +202,30 @@ END_SHELL
      enddo
 
      psi_value_save(i_walk) = psi_value
-     E_loc_save(i_walk) = E_loc
+     if (accepted) then
+        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
+     endif
 
 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_srmc_block_walk    += $X * srmc_pop_weight_mult * srmc_weight(i_walk)
-   !  $X_2_srmc_block_walk  += $X_2 * pop_weight_mult * srmc_weight(i_walk)
+   !  $X_2_srmc_block_walk  += $X_2 * srmc_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 3) = $X * srmc_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 3) = $X_2 * srmc_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)
@@ -203,7 +243,7 @@ for p in properties:
 
 END_SHELL
 
-    block_weight += pop_weight_mult * srmc_weight(i_walk)
+    block_weight += srmc_pop_weight_mult * srmc_weight(i_walk)
 
    else
      srmc_weight(i_walk) = 0.d0
@@ -221,15 +261,15 @@ END_SHELL
 
   ! Eventually, recompute the weight of the population
   if (srmc_projection_step == 1) then
-    pop_weight_mult = 1.d0
+    srmc_pop_weight_mult = 1.d0
     do k=1,srmc_projection
-      pop_weight_mult *= pop_weight(k)
+      srmc_pop_weight_mult *= srmc_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)
+  srmc_pop_weight_mult *= 1.d0/srmc_pop_weight(srmc_projection_step)
 
   ! Compute the new weight of the population
   double precision :: sum_weight
@@ -237,10 +277,10 @@ END_SHELL
   do k=1,walk_num
     sum_weight += srmc_weight(k)
   enddo
-  pop_weight(srmc_projection_step) = sum_weight/dble(walk_num)
+  srmc_pop_weight(srmc_projection_step) = sum_weight/dble(walk_num)
 
   ! Update the running population weight
-  pop_weight_mult *= pop_weight(srmc_projection_step)
+  srmc_pop_weight_mult *= srmc_pop_weight(srmc_projection_step)
 
 ! Reconfiguration
   integer :: ipos(walk_num)
@@ -257,7 +297,7 @@ END_SHELL
      enddo
     enddo
     psi_value_save_tmp(k) = psi_value_save(k)
-    E_loc_save_tmp(k) = E_loc_save(k)
+    E_loc_save_tmp(:,k) = E_loc_save(:,k)
   enddo
   
   integer :: ipm
@@ -272,7 +312,7 @@ END_SHELL
     enddo
    enddo
    psi_value_save(k) = psi_value_save_tmp(ipm)
-   E_loc_save(k) = E_loc_save_tmp(ipm)
+   E_loc_save(:,k) = E_loc_save_tmp(:,ipm)
   enddo
 
   call system_clock(cpu1, count_rate, count_max)
@@ -287,7 +327,7 @@ END_SHELL
     cpu2 = cpu1
   endif
 
-  SOFT_TOUCH elec_coord_full pop_weight_mult
+  SOFT_TOUCH elec_coord_full srmc_pop_weight_mult
 
   first_loop = .False.
 
@@ -314,7 +354,7 @@ END_SHELL
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, pop_weight_mult ]
+BEGIN_PROVIDER [ double precision, srmc_pop_weight_mult ]
  implicit none
  BEGIN_DOC  
 ! Population weight of SRMC
@@ -335,12 +375,12 @@ END_PROVIDER
  srmc_projection_step = 0
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, pop_weight, (0:srmc_projection+1) ]
+BEGIN_PROVIDER [ double precision, srmc_pop_weight, (0:srmc_projection+1) ]
  implicit none
  BEGIN_DOC  
 ! Population weight of SRMC
  END_DOC
- pop_weight = 1.d0
+ srmc_pop_weight = 1.d0
 END_PROVIDER