Merge pull request #2 from scemama/develop

Merge develop branch

ezfio_config/qmc.config

@@ -33,7 +33,6 @@ nuclei
   nucl_label          character*(32) (nuclei_nucl_num)
   nucl_charge         real           (nuclei_nucl_num)
   nucl_coord          real           (nuclei_nucl_num,3)
-  nucl_fitcusp_radius real           (nuclei_nucl_num)
 
 spindeterminants
   n_det_alpha     integer
@@ -55,7 +54,7 @@ simulation
   equilibration       logical
   http_server         character*(128)
   do_jast             logical
-  do_nucl_fitcusp     logical
+  nucl_fitcusp_factor real           
   method              character*(32)
   block_time          integer
   sampling            character*(32)

install/scripts/install_ocaml.sh

@@ -4,7 +4,7 @@ set -u
 set -e
 
 cd .. ; QMCCHEM_PATH="$PWD" ; cd -
-PACKAGES="core cryptokit ocamlfind sexplib"
+PACKAGES="core cryptokit ocamlfind sexplib pa_sexp_conv"
 
 declare -i i
 i=$(gcc -dumpversion | cut -d '.' -f 2)

ocaml/Default.ml

@@ -1,15 +1,21 @@
 open Core.Std;;
 
 
-let simulation_do_nucl_fitcusp  = lazy(
+let simulation_nucl_fitcusp_factor  = lazy(
-  if (not (Ezfio.has_pseudo_do_pseudo ())) then
+  let default = 
-    not (Ezfio.get_pseudo_do_pseudo ())
+    1.
+  in
+  if (Ezfio.has_pseudo_do_pseudo ()) then
+    if (Ezfio.get_pseudo_do_pseudo ()) then
+      0.
     else
-    true
+      default
+  else
+    default
 )
 
-let electrons_elec_walk_num         = lazy ( 30          )
+let electrons_elec_walk_num         = lazy ( 100         )
-let electrons_elec_walk_num_tot     = lazy ( 10000       )
+let electrons_elec_walk_num_tot     = lazy ( 1000        )
 let jastrow_jast_type               = lazy ( "None"      )
 let simulation_block_time           = lazy ( 30          )
 let simulation_ci_threshold         = lazy ( 1.e-8       )
@@ -26,9 +32,9 @@ let reset_defaults () =
       "/jastrow/jast_type"             ;
       "/simulation/block_time"         ;
       "/simulation/ci_threshold"       ;
-      "/simulation/do_nucl_fitcusp"    ;
       "/simulation/method"             ;
       "/simulation/sampling"           ;
       "/simulation/stop_time"          ;
-      "/simulation/time_step"          ]
+      "/simulation/time_step"          ;
+      "/simulation/nucl_fitcusp_factor" ]
 

ocaml/Input.ml

@@ -66,81 +66,69 @@ end = struct
 
 end
 
-module Fitcusp : sig
 
-  type t = bool
+
+module Fitcusp_factor : sig
+
+  type t = float
   val doc   : string
   val read  : unit -> t
   val write : t -> unit
-  val to_bool : t -> bool
+  val to_float : t -> float
-  val of_bool : bool -> t
+  val of_float : float -> t
-  val to_int  : t -> int
-  val of_int  : int -> t
   val to_string : t -> string
   val of_string : string -> t
 
 end = struct
 
-  type t = bool
+  type t = float
 
-  let doc = "Correct wave function to verify electron-nucleus cusp condition"
+  let doc = "Correct wave function to verify electron-nucleus cusp condition. 
+Fit is done for r < r_c(f) where r_c(f) = (1s orbital radius) x f. Value of f"
 
-  let of_bool x = x 
+  let of_float x = 
+    if (x < 0.) then
+       failwith "Fitcusp_factor should be >= 0.";
+    if (x > 10.) then
+       failwith "Fitcusp_factor is too large.";
+    x 
 
-  let to_bool x = x 
+  let to_float x = x 
 
   let read () = 
-    let _ =
+    ignore @@
-      Lazy.force Qputils.ezfio_filename
+      Lazy.force Qputils.ezfio_filename ;
+    if (not (Ezfio.has_simulation_nucl_fitcusp_factor ())) then
+      begin
+        let factor =
+          Lazy.force Default.simulation_nucl_fitcusp_factor ;
         in
-    if (not (Ezfio.has_simulation_do_nucl_fitcusp ())) then
+        Ezfio.set_simulation_nucl_fitcusp_factor factor
-      Lazy.force Default.simulation_do_nucl_fitcusp
+      end ;
-      |> Ezfio.set_simulation_do_nucl_fitcusp ;
+    Ezfio.get_simulation_nucl_fitcusp_factor () 
-    Ezfio.get_simulation_do_nucl_fitcusp ()
+    |> of_float 
-    |> of_bool
 
 
   let write t =
     let _ =
       Lazy.force Qputils.ezfio_filename
     in
-    let () =
+    to_float t
-      match (Pseudo.read () |> Pseudo.to_bool, to_bool t) with
+    |> Ezfio.set_simulation_nucl_fitcusp_factor 
-      | (true, true) -> failwith "Pseudopotentials and Fitcusp are incompatible"
-      | _ -> ()
-    in
-    to_bool t
-    |> Ezfio.set_simulation_do_nucl_fitcusp
 
 
   let to_string t =
-    to_bool t
+    to_float t
-    |> Bool.to_string
+    |> Float.to_string
 
 
   let of_string t =
     try
-      String.lowercase t
+      Float.of_string t
-      |> Bool.of_string 
+      |> of_float
-      |> of_bool
     with
     | Invalid_argument msg -> failwith msg
 
-
-  let to_int t =
-    let t = 
-      to_bool t
-    in
-    if t then 1
-    else 0
-
-
-  let of_int = function
-    | 0 -> false
-    | 1 -> true
-    | _ -> failwith "Expected 0 or 1"
-
-
 end
 
 module Block_time : sig
@@ -399,7 +387,7 @@ end
 
 module Method : sig
 
-  type t = VMC | DMC | SRMC
+  type t = VMC | DMC | SRMC | FKMC
   val doc : string
   val read  : unit -> t
   val write : t -> unit
@@ -408,21 +396,23 @@ module Method : sig
 
 end = struct
 
-  type t = VMC | DMC | SRMC
+  type t = VMC | DMC | SRMC | FKMC
 
-  let doc = "QMC Method : [ VMC | DMC | SRMC ]"
+  let doc = "QMC Method : [ VMC | DMC | SRMC | FKMC ]"
 
   let of_string = function
   | "VMC"  | "vmc"  -> VMC
   | "DMC"  | "dmc"  -> DMC
   | "SRMC" | "srmc" -> SRMC
-  | x -> failwith ("Method should be [ VMC | DMC | SRMC ], not "^x^".")
+  | "FKMC" | "fkmc" -> FKMC
+  | x -> failwith ("Method should be [ VMC | DMC | SRMC | FKMC ], not "^x^".")
 
 
   let to_string = function
   | VMC  -> "VMC"
   | DMC  -> "DMC"
   | SRMC -> "SRMC"
+  | FKMC -> "FKMC"
 
 
   let read () = 
@@ -855,8 +845,6 @@ let validate () =
      Time_step.read ()
   and jast_type = 
     Jastrow_type.read ()
-  and do_fitcusp =
-    Fitcusp.read ()
   and do_pseudo =
     Pseudo.read ()
   in
@@ -865,10 +853,12 @@ let validate () =
   let () =
     match (sampling, meth, Pseudo.to_bool do_pseudo) with
     | (Sampling.Brownian, Method.DMC, true) 
+    | (Sampling.Brownian, Method.FKMC, true) 
     | (Sampling.Brownian, Method.SRMC, true) ->
       if ( (Time_step.to_float ts) >= 0.5 ) then
           warn ( "Time step seems large for "^(Method.to_string meth) )
     | (Sampling.Brownian, Method.SRMC, false) 
+    | (Sampling.Brownian, Method.FKMC, false) 
     | (Sampling.Brownian, Method.DMC, false) ->
       if ( (Time_step.to_float ts) >= 0.01 ) then
           warn ( "Time step seems large for "^(Method.to_string meth) )
@@ -879,8 +869,9 @@ let validate () =
       if ( (Time_step.to_float ts) <= 0.01 ) then
           warn "Time step seems small for Langevin sampling."
     | (Sampling.Langevin, Method.SRMC, _) 
+    | (Sampling.Langevin, Method.FKMC, _) 
     | (Sampling.Langevin, Method.DMC, _) ->
-        failwith "Lanvegin sampling is incompatible with DMC/SRMC"
+        failwith "Lanvegin sampling is incompatible with DMC"
   in
 
 
@@ -888,6 +879,7 @@ let validate () =
   let () =
     match (meth, Ref_energy.(read () |> to_float) ) with
     | (Method.SRMC,0.) 
+    | (Method.FKMC,0.) 
     | (Method.DMC,0.) -> failwith ("E_ref should not be zero in "^(Method.to_string meth) )
     | _          -> ()
   in
@@ -902,6 +894,7 @@ let validate () =
   let () =
     match (meth, Property.(calc E_loc)) with
     | (Method.SRMC, false) 
+    | (Method.FKMC, false) 
     | (Method.DMC, false) -> failwith ( "E_loc should be sampled in "^(Method.to_string meth) )
     | (Method.VMC, false) -> warn "Sampling of E_loc is not activated in input"
     | _ -> ()
@@ -915,13 +908,23 @@ let validate () =
     | _ -> ()
   in
 
-  (* Fitcusp is not recommended with pseudo *)
+  (* Fitcusp is incompatible with pseudo *)
   let () =
-    match (Pseudo.to_bool do_pseudo, Fitcusp.to_bool do_fitcusp) with
+    let f = 
-    | (true, true) -> warn "Fitcusp is incompatible with Pseudopotentials"
+      Fitcusp_factor.read ()
+      |> Fitcusp_factor.to_float
+    in
+    match (Pseudo.to_bool do_pseudo, f > 0.) with
+    | (true, true) ->
+      begin
+        warn "Electron-nucleus cusp fitting is incompatible with Pseudopotentials.";
+        Fitcusp_factor.of_float 0.
+        |> Fitcusp_factor.write
+      end
     | _ -> ()
   in
 
+  
   (* Other Checks *)
   let () =
     let _ =

ocaml/Md5.ml

@@ -37,10 +37,9 @@ let files_to_track = [
   "mo_basis/mo_tot_num" ;
   "nuclei/nucl_charge.gz" ;
   "nuclei/nucl_coord.gz" ;
-  "nuclei/nucl_fitcusp_radius.gz" ;
   "nuclei/nucl_num" ;
   "simulation/ci_threshold" ;
-  "simulation/do_nucl_fitcusp" ;
+  "simulation/nucl_fitcusp_factor" ;
   "simulation/jast_a_up_dn" ;
   "simulation/jast_a_up_up" ;
   "simulation/jast_b_up_dn" ;

ocaml/Qmcchem_dataserver.ml

@@ -71,8 +71,6 @@ let run ?(daemon=true) ezfio_filename =
 
 
 
-  (*
-  (** Checks if the port is already open (not working properly yet) *)
   let check_port n =
     let adress_prefix =
       "tcp://*:"
@@ -87,7 +85,9 @@ let run ?(daemon=true) ezfio_filename =
         in
         let result = 
           try
-            (ZMQ.Socket.bind socket address; accu );
+            ZMQ.Socket.bind socket address;
+            ZMQ.Socket.unbind socket address; 
+            accu;
           with
           | _ -> false;
         in
@@ -100,22 +100,18 @@ let run ?(daemon=true) ezfio_filename =
     else
         `Unavailable
   in
-  *)
 
 
 
   (** Random port number between 49152 and 65535 *)
   let port =
     let newport = 
-  (*    ref (49152 + (Random.int 16383)) *)
       ref ( 1024 + (Random.int (49151-1024)))
     in
-  (*
     while ((check_port !newport) = `Unavailable)
     do
-      newport := 49152 + (Random.int 16383)
+      newport := 1024 + (Random.int (49151-1024))
     done;
-  *)
     !newport
   in
 
@@ -837,12 +833,13 @@ let run ?(daemon=true) ezfio_filename =
 
   (* Handle signals *)
   let handler s = 
+    Printf.printf "Dataserver received the %s signal... killing\n%!" (Signal.to_string s);
     Watchdog.kill ();
   in
   List.iter [
+     Signal.int  ;
      Signal.term ;
      Signal.quit ;
-     Signal.int  
     ]
     ~f:(fun x -> Signal.Expert.handle x handler)
   ;

ocaml/Qmcchem_edit.ml

@@ -20,7 +20,7 @@ type field =
  | Walk_num
  | Walk_num_tot
  | Stop_time
- | Fitcusp
+ | Fitcusp_factor
  | Method
  | Sampling
  | Ref_energy
@@ -54,8 +54,8 @@ let get field =
    option_to_string Walk_num_tot.read Walk_num_tot.to_string Walk_num_tot.doc  
  | Stop_time    ->
    option_to_string Stop_time.read    Stop_time.to_string    Stop_time.doc  
- | Fitcusp      ->
+ | Fitcusp_factor ->
-   option_to_string Fitcusp.read      Fitcusp.to_string      Fitcusp.doc
+   option_to_string Fitcusp_factor.read      Fitcusp_factor.to_string      Fitcusp_factor.doc
  | Method       ->
    option_to_string Method.read       Method.to_string       Method.doc  
  | Sampling     ->
@@ -136,7 +136,7 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename =
   handle_option         Input.Method.(of_string, write) m;
   handle_option      Input.Stop_time.(of_int   , write) t;
   handle_option       Input.Sampling.(of_string, write) s;
-  handle_option      Input.Fitcusp.(of_int   , write) f;
+  handle_option Input.Fitcusp_factor.(of_float , write) f;
   handle_option      Input.Time_step.(of_float , write) ts;
   handle_option       Input.Walk_num.(of_int   , write) w;
   handle_option   Input.Walk_num_tot.(of_int   , write) wt;
@@ -155,7 +155,7 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename =
       Ref_energy           ;
       Walk_num             ;
       Walk_num_tot         ;
-      Fitcusp              ;
+      Fitcusp_factor       ;
       CI_threshold         ;
       Jastrow_type         ;
       Properties           ;
@@ -214,7 +214,7 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename =
          begin
             match f with
             |  Stop_time             ->  Stop_time.(of_string            s  |>  write)
-            |  Fitcusp               ->  Fitcusp.(of_string              s  |>  write)
+            |  Fitcusp_factor        ->  Fitcusp_factor.(of_string       s  |>  write)
             |  Block_time            ->  Block_time.(of_string           s  |>  write)
             |  Method                ->  Method.(of_string               s  |>  write)
             |  Ref_energy            ->  Ref_energy.(of_string           s  |>  write)
@@ -271,8 +271,8 @@ let spec =
   empty
   +> flag "c"  no_arg
      ~doc:(" Clear blocks")
-  +> flag "f"  (optional int)
+  +> flag "f"  (optional float)
-     ~doc:("0|1 "^Input.Fitcusp.doc)
+     ~doc:("float "^Input.Fitcusp_factor.doc)
   +> flag "t"  (optional int)
      ~doc:("seconds "^Input.Stop_time.doc)
   +> flag "l"  (optional int)

ocaml/Qmcchem_forwarder.ml

@@ -1,14 +1,20 @@
 open Core.Std;;
 
 let bind_socket ~socket_type ~socket ~address =
+  let rec loop = function 
+  | 0 -> failwith @@ Printf.sprintf 
+        "Unable to bind the forwarder's %s socket : %s\n"
+        socket_type address
+  | -1 -> () 
+  | i ->  
       try 
-    ZMQ.Socket.bind socket address
+        ZMQ.Socket.bind socket address; 
+        loop (-1) 
       with 
-  | Unix.Unix_error (_, message, f) -> 
+      | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_float 1. ; loop (i-1) ) 
-    failwith @@ Printf.sprintf 
-        "\n%s\nUnable to bind the forwarder's %s socket :\n %s\n%s"
-        f socket_type address message
       | other_exception -> raise other_exception 
+  in loop 10 
+
 
 
 let run ezfio_filename dataserver =
@@ -45,18 +51,41 @@ let run ezfio_filename dataserver =
    *)
   let () = 
     try
-      Unix.mkdir tmpdir
+      Unix.mkdir tmpdir;
+      Unix.chdir tmpdir
     with
     | Unix.Unix_error _ ->
-        (* TODO : wait until the forwarder has started *)
         begin
           Unix.chdir tmpdir;
-          ignore @@ Unix.exec ~prog ~args ()
+          Time.pause @@ Time.Span.of_float 0.1;
+          match (Sys.file_exists "PID") with
+          | `No 
+          | `Unknown -> ()
+          | `Yes ->
+              let pid = 
+                In_channel.with_file "PID" ~f:(fun ic ->
+                  match (In_channel.input_line ic) with
+                  | Some x -> x
+                  | None -> "-1" )
+                |> Int.of_string
+              in
+              match pid with
+              | -1 -> ()
+              | pid ->
+                  begin
+                    match Signal.send (Signal.of_system_int 0) (`Pid (Pid.of_int pid)) with
+                    | `No_such_process ->  ()
+                    | _ -> ignore @@ Unix.exec ~prog ~args ()
+                  end
         end
   in
-  Unix.chdir tmpdir;
 
   (* Now, only one forwarder will execute the following code *)
+  Out_channel.with_file "PID" ~f:(fun oc ->
+    Unix.getpid ()
+    |> Pid.to_int
+    |> Printf.sprintf "%d\n" 
+    |> Out_channel.output_string oc);
 
   (* Fork a qmc *)
   ignore @@
@@ -90,20 +119,20 @@ let run ezfio_filename dataserver =
       with
       | _ ->  ()
       ;
-    done
+    done;
+    Watchdog.kill ()
   in
 
 
   (* Signal handler to Kill properly all the processes *)
   let handler s =
-    Printf.printf "Forwarder received the %s signal... killing\n" (Signal.to_string s);
+    Printf.printf "Forwarder received the %s signal... killing\n%!" (Signal.to_string s);
     terminate ();
-    Watchdog.kill ();
   in
   List.iter [
+     Signal.int  ;
      Signal.term ;
      Signal.quit ;
-     Signal.int
     ]
     ~f:(fun x -> Signal.Expert.handle x handler)
   ;
@@ -465,6 +494,13 @@ let run ezfio_filename dataserver =
   end;
 
   (* Wait for the qmc process to complete *)
+  try
     ignore (Watchdog.join ());
     terminate ()
+  with
+  | error -> 
+    begin
+      terminate ();
+      raise error
+    end
 

ocaml/Qmcchem_run.ml

@@ -7,8 +7,8 @@ let full_run ?(start_dataserver=true) ezfio_filename  =
   and scheduler = 
     Scheduler.find ()
   in
-  Printf.printf "Scheduler : %s\n" (Scheduler.to_string scheduler);
+  Printf.printf "Scheduler : %s\n%!" (Scheduler.to_string scheduler);
-  Printf.printf "Launcher  : %s\n" (Launcher.to_string  launcher );
+  Printf.printf "Launcher  : %s\n%!" (Launcher.to_string  launcher );
 
 
   (* Create the node file *)
@@ -147,13 +147,13 @@ let run a d ?q ?s ezfio_filename =
 
   (* Signal handler to Kill properly all the processes *)
   let handler s = 
-    Printf.printf "Received the %s signal... killing\n" (Signal.to_string s);
+    Printf.printf "QMC=Chem received the %s signal... killing\n%!" (Signal.to_string s);
     Watchdog.kill ();
   in
   List.iter [
+     Signal.int  ;
      Signal.term ;
      Signal.quit ;
-     Signal.int  
     ]
     ~f:(fun x -> Signal.Expert.handle x handler)
   ;

ocaml/Watchdog.ml

@@ -8,7 +8,7 @@ let _threads = ref [] ;;
 let kill () =
   let kill pid = 
     Signal.send_i Signal.int (`Pid pid);
-    Printf.printf "Killed %d\n" (Pid.to_int pid)
+    Printf.printf "Killed %d\n%!" (Pid.to_int pid)
   in
   List.iter ~f:kill (!_list);
   exit 1 

ocaml/build.ninja

@@ -1,7 +1,7 @@
 MAIN=qmcchem
 # Main program to build
 
-PACKAGES=-package core,sexplib.syntax,cryptokit,str,ZMQ
+PACKAGES=-package core,pa_sexp_conv,cryptokit,str,ZMQ
 # Required opam packages, for example:
 # PACKAGES=-package core,sexplib.syntax
 

src/AO/ao_full.irp.f

@@ -73,6 +73,7 @@ BEGIN_PROVIDER [ logical, primitives_reduced ]
       PROVIDE  ao_power
       PROVIDE  ao_coef
       PROVIDE  ao_nucl
+      PROVIDE  mo_fitcusp_normalization_before
       do i=1,ao_num
         if (ao_oned_p(i) /= 0.) then
           l=ao_power(i,1)+ao_power(i,2)+ao_power(i,3)

src/SAMPLING/block.irp.f

@@ -30,6 +30,12 @@ t = """
     $X_block_walk = $X_srmc_block_walk
     $X_2_block_walk = $X_2_srmc_block_walk
   endif
+ else if (qmc_method == t_FKMC) then
+  PROVIDE E_loc_fkmc_block_walk
+  if (calc_$X) then
+    $X_block_walk = $X_fkmc_block_walk
+    $X_2_block_walk = $X_2_fkmc_block_walk
+  endif
  endif
 
 END_PROVIDER

src/SAMPLING/dmc_step.irp.f

@@ -195,7 +195,7 @@ END_SHELL
   do k=1,walk_num_dmc
     sum_weight += dmc_weight(k)
   enddo
-  E0 = E_ref - log(real(walk_num_dmc)/real(walk_num)) * 0.1d0/dtime_step
+  E0 = E_ref - log(sum_weight/real(walk_num)) * 0.1d0 /dtime_step
 
 ! Branching
   integer                        :: ipos(walk_num_dmc_max), walk_num_dmc_new
@@ -214,20 +214,22 @@ END_SHELL
     ipos(k) = k
   enddo
 
+  walk_num_dmc_new = walk_num_dmc
   do k=1,walk_num_dmc
     r = qmc_ranf()
     if (dmc_weight(k) > 1.d0) then
       if ( 1.d0+r < dmc_weight(k) ) then
-        walk_num_dmc = walk_num_dmc+1
+        walk_num_dmc_new = walk_num_dmc_new+1
-        ipos(walk_num_dmc) = k
+        ipos(walk_num_dmc_new) = k
       endif
     else
       if ( r > dmc_weight(k) ) then
-        ipos(k) = ipos(walk_num_dmc)
+        ipos(k) = ipos(walk_num_dmc_new)
-        walk_num_dmc = walk_num_dmc-1
+        walk_num_dmc_new = walk_num_dmc_new-1
       endif
     endif
   enddo
+  walk_num_dmc = walk_num_dmc_new
 
   integer :: ipm
   do k=1,walk_num_dmc
@@ -328,7 +330,7 @@ BEGIN_PROVIDER [ integer, walk_num_dmc_max ]
  BEGIN_DOC
  ! Max number of walkers in DMC
  END_DOC
- walk_num_dmc_max = 3 * walk_num
+ walk_num_dmc_max = max(3 * walk_num, 30)
 END_PROVIDER
 
 

src/SAMPLING/fkmc_step.irp.f

@@ -0,0 +1,375 @@
+! Providers of *_fkmc_block_walk
+!==============================
+BEGIN_SHELL [ /usr/bin/python ]
+from properties import *
+
+t = """
+ BEGIN_PROVIDER [ $T, $X_fkmc_block_walk $D1 ]
+&BEGIN_PROVIDER [ $T, $X_fkmc_block_walk_kahan $D2 ]
+&BEGIN_PROVIDER [ $T, $X_2_fkmc_block_walk $D1 ]
+&BEGIN_PROVIDER [ $T, $X_2_fkmc_block_walk_kahan $D2 ]
+ implicit none
+ BEGIN_DOC  
+! fkMC averages of $X. Computed in E_loc_fkmc_block_walk
+ END_DOC
+ $X_fkmc_block_walk = 0.d0
+ $X_fkmc_block_walk_kahan = 0.d0
+ $X_2_fkmc_block_walk = 0.d0
+ $X_2_fkmc_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_fkmc_block_walk       ]
+&BEGIN_PROVIDER [ double precision, E_loc_2_fkmc_block_walk     ]
+&BEGIN_PROVIDER [ double precision, E_loc_fkmc_block_walk_kahan, (3) ]
+&BEGIN_PROVIDER [ double precision, E_loc_2_fkmc_block_walk_kahan, (3) ]
+ implicit none
+ include '../types.F'
+ BEGIN_DOC
+! Properties averaged over the block using the FKMC method
+ END_DOC
+
+  integer, parameter :: BIRTH=1, DEATH=2
+  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 :: fkmc_weight(walk_num)
+  double precision :: delta(walk_num)
+  double precision, allocatable :: psi_grad_psi_inv_save(:,:,:)
+  double precision, allocatable :: psi_grad_psi_inv_save_tmp(:,:,:)
+  double precision, allocatable :: fkmc_clock_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 :: 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 :: fkmc_weight
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: fkmc_clock_tmp
+  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),             &
+       fkmc_clock_tmp(2,walk_num) )
+  psi_value_save = 0.d0
+  psi_value_save_tmp = 0.d0
+  fkmc_weight = 1.d0
+
+! Initialization
+ if (vmc_algo /= t_Brownian) then
+   call abrt(irp_here,'FKMC 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_fkmc_block_walk = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_fkmc_block_walk_kahan = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_2_fkmc_block_walk = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_2_fkmc_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               :: thr
+
+ thr = 2.d0/time_step_sq
+
+ logical :: first_loop
+ first_loop = .True.
+
+ do while (loop)
+
+  ! Every walker makes a step
+  do i_walk=1,walk_num
+    
+    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(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
+      E_loc_save(i_walk) = E_loc 
+      psi_value_save(i_walk) = psi_value
+    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
+     delta(i_walk) = ((E_loc+E_loc_save(i_walk))*0.5d0 - E_ref) * p
+     if ( delta(i_walk) > thr ) then
+       delta(i_walk) = thr
+     else if ( delta(i_walk) < -thr ) then
+       delta(i_walk) = -thr
+     endif
+     fkmc_weight(i_walk) = dexp(-dtime_step*delta(i_walk))
+     elec_coord(elec_num+1,1) += p*time_step
+     elec_coord(elec_num+1,2)  = E_loc
+     elec_coord(elec_num+1,3)  = fkmc_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
+     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(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_fkmc_block_walk    += $X * fkmc_weight(i_walk)
+   !  $X_2_fkmc_block_walk  += $X_2 * fkmc_weight(i_walk)
+   ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
+   
+      $X_fkmc_block_walk_kahan($D2 3) = $X * fkmc_weight(i_walk) - $X_fkmc_block_walk_kahan($D2 1)
+      $X_fkmc_block_walk_kahan($D2 2) = $X_fkmc_block_walk $D1  + $X_fkmc_block_walk_kahan($D2 3)
+      $X_fkmc_block_walk_kahan($D2 1) = ($X_fkmc_block_walk_kahan($D2 2) - $X_fkmc_block_walk $D1 ) &
+          - $X_fkmc_block_walk_kahan($D2 3)
+      $X_fkmc_block_walk $D1  =  $X_fkmc_block_walk_kahan($D2 2) 
+   
+   
+      $X_2_fkmc_block_walk_kahan($D2 3) = $X_2 * fkmc_weight(i_walk) - $X_2_fkmc_block_walk_kahan($D2 1)
+      $X_2_fkmc_block_walk_kahan($D2 2) = $X_2_fkmc_block_walk $D1 + $X_2_fkmc_block_walk_kahan($D2 3)
+      $X_2_fkmc_block_walk_kahan($D2 1) = ($X_2_fkmc_block_walk_kahan($D2 2) - $X_2_fkmc_block_walk $D1 ) &
+          - $X_2_fkmc_block_walk_kahan($D2 3)
+      $X_2_fkmc_block_walk $D1 =  $X_2_fkmc_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 += fkmc_weight(i_walk)
+
+   else
+     fkmc_weight(i_walk) = 0.d0
+     delta(i_walk) = 1.d5
+   endif
+
+  enddo
+
+  ! Compute the new weight of the population
+  double precision :: sum_weight
+  sum_weight = 0.d0
+  do k=1,walk_num
+    sum_weight += fkmc_weight(k)
+  enddo
+
+  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
+     do i=1,elec_num
+      psi_grad_psi_inv_save_tmp(i,l,k) = psi_grad_psi_inv_save(i,l,k)
+     enddo
+    enddo
+    psi_value_save_tmp(k) = psi_value_save(k)
+    E_loc_save_tmp(k) = E_loc_save(k)
+    if (fkmc_weight(k) == 0.d0) then
+      fkmc_clock(DEATH,k) = -1.d0
+    endif
+    if ( delta(k) <= 0.d0 ) then
+      fkmc_clock_tmp(BIRTH,k) = fkmc_clock(BIRTH,k) +time_step * delta(k)
+      fkmc_clock_tmp(DEATH,k) = fkmc_clock(DEATH,k)
+    else
+      fkmc_clock_tmp(BIRTH,k) = fkmc_clock(BIRTH,k)
+      fkmc_clock_tmp(DEATH,k) = fkmc_clock(DEATH,k) -time_step * delta(k)
+    endif
+  enddo
+
+! Reconfiguration
+! ===============
+
+  ! Identify first which walkers will be killed to place branched walkers there
+  ! later
+  
+  double precision, external     :: qmc_ranf
+  integer                        :: ipm, m
+  integer                        :: killed(walk_num)
+
+  m=1
+  do k=1,walk_num
+    fkmc_clock(DEATH,k) = fkmc_clock_tmp(DEATH,k) 
+    if (fkmc_clock_tmp(DEATH,k) <= 0.d0) then
+      killed(m) = k
+      m += 1
+      fkmc_clock(DEATH,k) = -dlog(qmc_ranf())
+      fkmc_clock(BIRTH,k) = -dlog(qmc_ranf())
+      ipm = k
+      do while (ipm == k)
+        ipm = 1 + int (walk_num*qmc_ranf())
+      enddo
+      do l=1,3
+        do i=1,elec_num+1
+        elec_coord_full(i,l,k) = elec_coord_tmp(i,l,ipm)
+        enddo
+        do i=1,elec_num
+          psi_grad_psi_inv_save(i,l,k) = psi_grad_psi_inv_save_tmp(i,l,ipm)
+        enddo
+      enddo
+      psi_value_save(k) = psi_value_save_tmp(ipm)
+      E_loc_save(k) = E_loc_save_tmp(ipm)
+    endif
+  enddo
+  killed(m) = 0
+
+  m=1
+  do k=1,walk_num
+    fkmc_clock(BIRTH,k) = fkmc_clock_tmp(BIRTH,k) 
+    if (fkmc_clock_tmp(BIRTH,k) <= 0.d0) then
+      fkmc_clock(BIRTH,k) = -dlog(qmc_ranf())
+      if (killed(m) == 0) then
+        ipm = k
+        do while (ipm == k)
+          ipm = 1 + int (walk_num*qmc_ranf())
+        enddo
+      else
+        ipm = killed(m)
+        m +=1
+      endif
+      fkmc_clock(BIRTH,ipm) = -dlog(qmc_ranf())
+      fkmc_clock(DEATH,ipm) = -dlog(qmc_ranf())
+      do l=1,3
+        do i=1,elec_num+1
+        elec_coord_full(i,l,ipm) = elec_coord_tmp(i,l,k)
+        enddo
+        do i=1,elec_num
+          psi_grad_psi_inv_save(i,l,ipm) = psi_grad_psi_inv_save_tmp(i,l,k)
+        enddo
+      enddo
+      psi_value_save(ipm) = psi_value_save_tmp(k)
+      E_loc_save(ipm) = E_loc_save_tmp(k)
+      
+    endif
+
+  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
+
+! Update E_ref to take into account the weight of the population
+  E_ref -= dlog(sum_weight / dble(walk_num) ) / time_step 
+  SOFT_TOUCH elec_coord_full E_ref
+
+  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_fkmc_block_walk   *= factor
+   $X_2_fkmc_block_walk *= factor
+ endif
+"""
+for p in properties:
+ print  t.replace("$X",p[1])
+END_SHELL
+
+ deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp, &
+   fkmc_clock_tmp )
+
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [ double precision, fkmc_clock, (2,walk_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Branching clocks for the FKMC algotithm. (1,:) is the birth clock and
+ ! (2,:) is the death clock.
+ END_DOC
+ integer :: i
+ double precision, external :: qmc_ranf
+ do i=1, walk_num
+   fkmc_clock(1,i) = -dlog(qmc_ranf())
+   fkmc_clock(2,i) = -dlog(qmc_ranf())
+ enddo
+
+END_PROVIDER
+

src/SAMPLING/srmc_step.irp.f

@@ -133,7 +133,9 @@ END_SHELL
           elec_coord(i,l) = elec_coord_full(i,l,i_walk)
         enddo
       enddo
-      SOFT_TOUCH elec_coord
+      TOUCH elec_coord
+      psi_value_save(i_walk) = psi_value
+      E_loc_save(i_walk) = E_loc
     endif
 
    double precision               :: p,q

src/ezfio_interface.irp.f

@@ -6,7 +6,6 @@ data = [ \
 ("nuclei_nucl_num"                 , "integer"       , ""                      ),
 ("nuclei_nucl_charge"              , "real"          , "(nucl_num)"            ),
 ("nuclei_nucl_coord"               , "real"          , "(nucl_num,3)"          ),
-("nuclei_nucl_fitcusp_radius"      , "real"          , "(nucl_num)"            ),
 ("mo_basis_mo_coef"                , "real"          , "(ao_num,mo_tot_num)" ),
 ("electrons_elec_fitcusp_radius"   , "real"          , ""                      ),
 ("electrons_elec_alpha_num"        , "integer"       , ""                      ),
@@ -38,9 +37,9 @@ data = [ \
 ("simulation_time_step"            , "real"          , ""                      ),
 ("simulation_srmc_projection_time" , "real"          , ""                      ),
 ("simulation_method"               , "character*(32)", ""                      ),
+("simulation_nucl_fitcusp_factor"  , "real"          , ""                     ),
 ("simulation_save_data"            , "logical"       , ""                      ),
 ("simulation_print_level"          , "integer"       , ""                      ),
-("simulation_do_nucl_fitcusp"      , "logical"       , ""                      ),
 ("simulation_sampling"             , "character*(32)", ""                      ),
 ("simulation_ci_threshold"         , "double precision"          , ""                      ),
 ("simulation_http_server"          , "character*(128)", ""                      ),

src/mo.irp.f

@@ -273,6 +273,15 @@ END_PROVIDER
     enddo
   endif
   
+  do i=1,mo_num
+    do j=1,elec_num
+        mo_value_transp(i,j)  *= mo_cusp_rescale(i)
+        mo_grad_transp_x(i,j) *= mo_cusp_rescale(i)
+        mo_grad_transp_y(i,j) *= mo_cusp_rescale(i)
+        mo_grad_transp_z(i,j) *= mo_cusp_rescale(i)
+        mo_lapl_transp(i,j)   *= mo_cusp_rescale(i)
+    enddo
+  enddo
 END_PROVIDER
 
 
@@ -401,6 +410,7 @@ BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ]
   integer                        :: i, j, k, l
   real                           :: ao_value_at_nucl_no_S(ao_num)
   
+  PROVIDE mo_fitcusp_normalization_before
   do k=1,nucl_num
     point(1) = nucl_coord(k,1)
     point(2) = nucl_coord(k,2)
@@ -466,6 +476,99 @@ END_PROVIDER
   FREE ao_value_p ao_grad_p ao_lapl_p ao_axis_grad_p ao_oned_grad_p ao_oned_prim_grad_p ao_oned_lapl_p ao_axis_lapl_p ao_oned_prim_lapl_p ao_oned_p ao_oned_prim_p ao_axis_p ao_axis_power_p
   SOFT_TOUCH point
   
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_before, (mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Renormalization factor of MOs due to cusp fitting
+ END_DOC
+ include 'constants.F'
+ integer :: i,j,k,l
+ double precision :: dr, r, f, t
+ integer, save :: ifirst = 0
+
+  if (ifirst == 0) then
+    ifirst = 1
+    mo_fitcusp_normalization_before = 0.d0
+    do k=1,nucl_num
+      dr = nucl_fitcusp_radius(k)*1.d-2
+      point(1) = nucl_coord(k,1)
+      point(2) = nucl_coord(k,2)
+      point(3) = nucl_coord(k,3)-dr
+      do l=1,101
+        r = point(3) + dr
+        point(3) = r
+        TOUCH point
+        f = dfour_pi*r*r*dr
+        do i=1,mo_tot_num
+          mo_fitcusp_normalization_before(i) += f*mo_value_p(i)**2
+        enddo
+      enddo
+    enddo
+  endif
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_after, (mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Renormalization factor of MOs due to cusp fitting
+ END_DOC
+ include 'constants.F'
+ integer :: i,j,k,l
+ double precision :: dr, r, f, t, t2
+ integer, save :: ifirst = 0
+
+  PROVIDE primitives_reduced
+  if (ifirst == 0) then
+    ifirst = 1
+    mo_fitcusp_normalization_after = 0.d0
+    do k=1,nucl_num
+      dr = nucl_fitcusp_radius(k)*1.d-2
+      point(1) = nucl_coord(k,1)
+      point(2) = nucl_coord(k,2)
+      point(3) = nucl_coord(k,3)- dr
+      do l=1,101
+        point(3) = point(3)+ dr
+        TOUCH point nucl_fitcusp_param primitives_reduced mo_coef
+        r = point(3)
+        f = dfour_pi*r*r*dr
+        do i=1,mo_num
+          t = 0.d0
+          do j=1,ao_num
+            if ( (ao_nucl(j) /= k).or.(ao_power(j,4) > 0) ) then
+              t = t + mo_coef(j,i) * ao_value_p(j) 
+            endif
+          enddo
+          t = t + nucl_fitcusp_param(1,i,k) +                        &
+              r * (nucl_fitcusp_param(2,i,k) +                       &
+              r * (nucl_fitcusp_param(3,i,k) +                       &
+              r *  nucl_fitcusp_param(4,i,k) ))
+          mo_fitcusp_normalization_after(i) += t*t*f
+        enddo
+      enddo
+    enddo
+  endif
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ real, mo_cusp_rescale, (mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Rescaling coefficient to normalize MOs after applying fitcusp
+ END_DOC
+ integer :: i
+ if (do_nucl_fitcusp) then
+   do i=1,mo_tot_num
+!     mo_cusp_rescale(i) = dsqrt(mo_fitcusp_normalization_before(i) / mo_fitcusp_normalization_after(i))
+     mo_cusp_rescale(i) = 1.d0/dsqrt(1.d0 - mo_fitcusp_normalization_before(i) + mo_fitcusp_normalization_after(i))
+   enddo
+ else
+     mo_cusp_rescale = 1.d0
+ endif
+
 END_PROVIDER
 
 

src/nuclei.irp.f

@@ -125,17 +125,20 @@ BEGIN_PROVIDER [ real, nucl_fitcusp_radius, (nucl_num) ]
   BEGIN_DOC
 ! Distance threshold for the fit
   END_DOC
-  real                           :: def(nucl_num)
+  real                           :: def(nucl_num), factor
   integer                        :: k
+  real, parameter                :: a = 1.74891
+  real, parameter                :: b = 0.126057
   
-  if (.not.do_nucl_fitcusp) then
+  if (.not. do_nucl_fitcusp) then
     nucl_fitcusp_radius = 0.d0
     return
   endif
+
   do k=1,nucl_num
-    nucl_fitcusp_radius(k) = .5/nucl_charge(k)
+    nucl_fitcusp_radius(k) = nucl_fitcusp_factor/(a*nucl_charge(k)+b)
   enddo
-  call get_nuclei_nucl_fitcusp_radius(nucl_fitcusp_radius)
+
   ! Avoid dummy atoms
   do k=1,nucl_num
     if (nucl_charge(k) < 5.d-1) then

src/simulation.irp.f

@@ -148,7 +148,7 @@ BEGIN_PROVIDER  [ integer, qmc_method ]
    implicit none
    include 'types.F'
    BEGIN_DOC
-   ! qmc_method : Calculation method. Can be t_VMC, t_DMC, t_SRMC
+   ! qmc_method : Calculation method. Can be t_VMC, t_DMC, t_SRMC, t_FKMC
    END_DOC
    character*(32)                 :: method
    method = types(t_VMC)
@@ -160,8 +160,10 @@ BEGIN_PROVIDER  [ integer, qmc_method ]
      qmc_method = t_DMC
    else if (method == types(t_SRMC)) then
      qmc_method = t_SRMC
+   else if (method == types(t_FKMC)) then
+     qmc_method = t_FKMC
    else
-     call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC )')
+     call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC | FKMC )')
    endif
    
    call cinfo(irp_here,'qmc_method',trim(method))
@@ -250,17 +252,22 @@ BEGIN_PROVIDER [ character*(64), hostname]
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ logical, do_nucl_fitcusp ]
+ BEGIN_PROVIDER [ real, nucl_fitcusp_factor ]
+&BEGIN_PROVIDER [ logical, do_nucl_fitcusp ]
  implicit none
  BEGIN_DOC
-   ! If true, do the fit of the electron-nucleus cusp
+ ! The electron-nucleus cusp fitting is done between 0 and r_c,
+ ! where r_c is chosen as nucl_fitcusp_factor * (radius_of_1s AO)
  END_DOC
-   do_nucl_fitcusp = .True.
+ nucl_fitcusp_factor = 0.
-   call get_simulation_do_nucl_fitcusp(do_nucl_fitcusp)
+ call get_simulation_nucl_fitcusp_factor(nucl_fitcusp_factor)
-   call linfo(irp_here,'do_nucl_fitcusp',do_nucl_fitcusp)
+ do_nucl_fitcusp = nucl_fitcusp_factor > 0.
+ call info(irp_here,'nucl_fitcusp_factor',nucl_fitcusp_factor)
+
 END_PROVIDER
 
 
+ 
 BEGIN_PROVIDER [ integer, vmc_algo ]
    implicit none
    include 'types.F'
@@ -283,6 +290,9 @@ BEGIN_PROVIDER [ integer, vmc_algo ]
      if (qmc_method == t_SRMC) then
        stop 'Langevin incompatible with SRMC'
      endif
+     if (qmc_method == t_FKMC) then
+       stop 'Langevin incompatible with FKMC'
+     endif
    else if (Sampling == types(t_MTM)) then
      vmc_algo = t_MTM
    else

src/types.F

@@ -6,6 +6,7 @@
       integer, parameter :: t_VMC      = 7
       integer, parameter :: t_DMC      = 8
       integer, parameter :: t_SRMC     = 9
+      integer, parameter :: t_FKMC     = 10
 
       integer, parameter :: t_Simple   = 11
       integer, parameter :: t_None     = 12
@@ -26,7 +27,7 @@
           'VMC            ',   &
           'DMC            ',   &
           'SRMC           ',   &
-          '               ',   &
+          'FKMC           ',   &
           'Simple         ',   &
           'None           ',   &
           '               ',   &