Repaired broken DMC

ocaml/Input.ml

@@ -630,7 +630,7 @@ end = struct
   let of_float x = 
     if (x >= 100.) then
       failwith "DMC Projection time should be < 100.";
-    if (x < 0.) then
+    if (x <= 0.) then
       failwith "DMC Projection time should be positive.";
     x
 
@@ -886,16 +886,6 @@ let validate () =
     | _          -> ()
   in
 
-(*
-  (* Check Projection time is greater than time step *)
-  let () =
-    match (meth, DMC_projection_time.(read () |> to_float) ) with
-    | (Method.DMC,p) -> 
-      if (p < ts) then failwith "E_ref should not be zero in DMC"
-    | _          -> ()
-  in
-*)
-
   (* Set block and total time*)
   let () =
     if ( (Block_time.read ()) > Stop_time.read ()) then

src/SAMPLING/dmc_step.irp.f

@@ -135,6 +135,7 @@ for p in properties:
 END_SHELL
 
 
+   ! Brownian step
    double precision               :: p,q
    real                           :: delta_x
    logical                        :: accepted

src/SAMPLING/reconfigure.irp.f

@@ -1,3 +1,52 @@
+subroutine reconfigure_old(ipos,w)
+  implicit none
+  integer, intent(inout)         :: ipos(walk_num)
+  double precision, intent(in)   :: w(walk_num)
+  
+  integer                        :: kptab(walk_num), kmtab(walk_num)
+  double precision               :: wp(walk_num), wm(walk_num)
+  double precision               :: tmp
+  
+  integer                        :: k, l
+  double precision               :: qmc_ranf, rand
+  integer                        :: ipos_tmp(walk_num*2)
+  
+  l=0
+  do k=1,walk_num
+    tmp = w(k)-1.d0
+    rand = qmc_ranf()
+    if (tmp >= 0.d0) then
+      l=l+1
+      ipos_tmp(l) = k
+      if (rand < tmp) then
+        l=l+1
+        ipos_tmp(l) = k
+      endif
+    else
+      if (rand > -tmp) then
+        l=l+1
+        ipos_tmp(l) = k
+      endif
+    endif
+  enddo
+  if (l>walk_num) then
+    do k=l,walk_num+1,-1
+      rand = qmc_ranf() * dble(walk_num)
+      ipos_tmp(int(rand)+1) = ipos_tmp(k) 
+    enddo
+  else
+    do k=l+1,walk_num
+      rand = qmc_ranf() * dble(k-1)
+      ipos_tmp(k) = ipos_tmp(int(rand)+1)
+    enddo
+  endif
+
+  do k=1,walk_num
+    ipos(k) = ipos_tmp(k)
+  enddo
+
+end
+
 subroutine reconfigure(ipos,w)
   implicit none
   integer, intent(inout)         :: ipos(*)