diff --git a/install/build.ninja b/install/build.ninja
index 351db7d..34ab4f7 100644
--- a/install/build.ninja
+++ b/install/build.ninja
@@ -7,15 +7,14 @@ URL_OPAM  ="https://raw.github.com/ocaml/opam/master/shell/opam_installer.sh"
 URL_IRPF90="https://github.com/scemama/irpf90/archive/v1.6.7.tar.gz"
 URL_EZFIO ="https://github.com/scemama/EZFIO/archive/v1.3.1.tar.gz"
 
-URL_ZMQ   ="http://download.zeromq.org/zeromq-4.0.7.tar.gz"
-#URL_ZMQ   ="http://download.zeromq.org/zeromq-4.1.3.tar.gz"
-URL_F77ZMQ="https://github.com/scemama/f77_zmq/archive/v4.1.3.tar.gz"
+URL_ZMQ   ="http://download.zeromq.org/zeromq-4.1.4.tar.gz"
+URL_F77ZMQ="https://github.com/scemama/f77_zmq/archive/4.1.4.tar.gz"
 
 # Rules
 #######
 
 rule download
-  command = [[ -e ${out} ]] || (wget --no-check-certificate ${url} -O ${out}.tmp -o /dev/null && mv ${out}.tmp ${out})
+  command = [ -e ${out} ] || (wget --no-check-certificate ${url} -O ${out}.tmp -o /dev/null && mv ${out}.tmp ${out})
   description = Downloading ${descr}
 
 rule install
diff --git a/install/scripts/install_ocaml.sh b/install/scripts/install_ocaml.sh
index 5718ed3..e6f0f0b 100755
--- a/install/scripts/install_ocaml.sh
+++ b/install/scripts/install_ocaml.sh
@@ -6,9 +6,43 @@ set -e
 cd .. ; QMCCHEM_PATH="$PWD" ; cd -
 PACKAGES="core cryptokit ocamlfind sexplib" # ppx_sexp_conv"
 
-declare -i i
-i=$(gcc -dumpversion | cut -d '.' -f 2)
-if [[ i -lt 6 ]]
+# return 0 if program version is equal or greater than check version
+check_version () {
+    if [[ $1 == $2 ]]
+    then
+        return 0
+    fi
+    local IFS=.
+    local i ver1=($1) ver2=($2)
+    # fill empty fields in ver1 with zeros
+    for ((i=${#ver1[@]}; i<${#ver2[@]}; i++))
+    do
+        ver1[i]=0
+    done
+    for ((i=0; i<${#ver1[@]}; i++))
+    do
+        if [[ -z ${ver2[i]} ]]
+        then
+            # fill empty fields in ver2 with zeros
+            ver2[i]=0
+        fi
+        if ((10#${ver1[i]} > 10#${ver2[i]}))
+        then
+            return 1
+        fi
+        if ((10#${ver1[i]} < 10#${ver2[i]}))
+        then
+            return 0
+        fi
+    done
+    return 0
+}
+
+
+i=$(gcc -dumpversion)
+
+check_version 4.6 $i
+if [[ $? == 1 ]]
 then
    echo "GCC version $(gcc -dumpversion) too old. GCC >= 4.6 required."
    exit 1
diff --git a/install/scripts/install_ocaml_zmq.sh b/install/scripts/install_ocaml_zmq.sh
index 1591fb0..a6f73d7 100755
--- a/install/scripts/install_ocaml_zmq.sh
+++ b/install/scripts/install_ocaml_zmq.sh
@@ -32,7 +32,7 @@ export C_INCLUDE_PATH="${QMCCHEM_PATH}/lib":$C_INCLUDE_PATH
 export LIBRARY_PATH="${QMCCHEM_PATH}/lib":$LIBRARY_PATH
 export LD_LIBRARY_PATH="${QMCCHEM_PATH}/lib":$LD_LIBRARY_PATH
 set -u
-opam install zmq 
+opam install zmq conf-zmq
 rm -f _build/ocaml_zmq.log
 exit 0
 
diff --git a/install/scripts/install_zmq.sh b/install/scripts/install_zmq.sh
index 80da3a3..6c7609c 100755
--- a/install/scripts/install_zmq.sh
+++ b/install/scripts/install_zmq.sh
@@ -3,7 +3,7 @@
 TARGET=zmq
 function _install()
 {
-  LIBVERSION=4
+  LIBVERSION=5
   cd .. ; QMCCHEM_PATH="$PWD" ; cd -
   set +u
   export C_INCLUDE_PATH="${C_INCLUDE_PATH}":./
@@ -14,10 +14,10 @@ function _install()
   make -j 8
   cd - 
   rm -f -- "${QMCCHEM_PATH}"/lib/libzmq.{a,so,so.$LIBVERSION} 
-# cp "${BUILD}"/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/ 
-# cp "${BUILD}"/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION	
-  cp "${BUILD}"/src/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/ 
-  cp "${BUILD}"/src/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION
+ cp "${BUILD}"/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/ 
+ cp "${BUILD}"/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION	
+#  cp "${BUILD}"/src/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/ 
+#  cp "${BUILD}"/src/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION
   cp "${BUILD}"/include/{zmq,zmq_utils}.h "${QMCCHEM_PATH}"/lib/
   cd "${QMCCHEM_PATH}"/lib
   ln libzmq.so.$LIBVERSION libzmq.so || cp libzmq.so.$LIBVERSION libzmq.so
diff --git a/ocaml/Block.ml b/ocaml/Block.ml
index b78eabe..a26c1f8 100644
--- a/ocaml/Block.ml
+++ b/ocaml/Block.ml
@@ -1,5 +1,5 @@
-open Core.Std;;
-open Qptypes;;
+open Core
+open Qptypes
 
 type t = 
   { property       : Property.t ;
diff --git a/ocaml/Default.ml b/ocaml/Default.ml
index 12f5f86..c0c3dce 100644
--- a/ocaml/Default.ml
+++ b/ocaml/Default.ml
@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core
 
 
 let simulation_nucl_fitcusp_factor  = lazy(
diff --git a/ocaml/Input.ml b/ocaml/Input.ml
index 70b321b..028acc2 100644
--- a/ocaml/Input.ml
+++ b/ocaml/Input.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 open Qputils
 
diff --git a/ocaml/Launcher.ml b/ocaml/Launcher.ml
index 4b9de5a..7b4066d 100644
--- a/ocaml/Launcher.ml
+++ b/ocaml/Launcher.ml
@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core
 
 type t =
   | Srun
diff --git a/ocaml/Md5.ml b/ocaml/Md5.ml
index 48d0b9b..fa177f1 100644
--- a/ocaml/Md5.ml
+++ b/ocaml/Md5.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 (** Directory containing the list of input files. The directory is created is inexistant. *)
 let input_directory = lazy (
diff --git a/ocaml/Message.ml b/ocaml/Message.ml
index 8e103ae..2db843b 100644
--- a/ocaml/Message.ml
+++ b/ocaml/Message.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 type t =
diff --git a/ocaml/Qmcchem_config.ml b/ocaml/Qmcchem_config.ml
index 894279a..d9f302d 100644
--- a/ocaml/Qmcchem_config.ml
+++ b/ocaml/Qmcchem_config.ml
@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core
 
 
 (** QMC=Chem installation directory *)
diff --git a/ocaml/Qmcchem_dataserver.ml b/ocaml/Qmcchem_dataserver.ml
index 3dce4f1..3067996 100644
--- a/ocaml/Qmcchem_dataserver.ml
+++ b/ocaml/Qmcchem_dataserver.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 (** Data server of QMC=Chem.
@@ -36,7 +36,7 @@ let run ?(daemon=true) ezfio_filename =
         begin
           Printf.printf "Generating initial walkers...\n%!";
           Unix.fork_exec ~prog:(Lazy.force Qmcchem_config.qmc_create_walkers)
-              ~args:["qmc_create_walkers" ; ezfio_filename] ()
+              ~argv:["qmc_create_walkers" ; ezfio_filename] ()
           |> Unix.waitpid_exn ;
           Printf.printf "Initial walkers ready\n%!"
         end ;
@@ -98,8 +98,7 @@ let run ?(daemon=true) ezfio_filename =
         let result = 
           try
             ZMQ.Socket.bind socket address;
-            ZMQ.Socket.unbind socket address; 
-            accu;
+            accu
           with
           | _ -> false;
         in
diff --git a/ocaml/Qmcchem_debug.ml b/ocaml/Qmcchem_debug.ml
index 950446e..887d90d 100644
--- a/ocaml/Qmcchem_debug.ml
+++ b/ocaml/Qmcchem_debug.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 
 let run ~t ezfio_filename=
@@ -49,7 +49,7 @@ let run ~t ezfio_filename=
           in
           tot_size := Byte_units.create `Bytes ((Byte_units.bytes !tot_size) +. (Byte_units.bytes bytes));
           Printf.printf "%s\n%!" (Byte_units.to_string !tot_size);
-          Time.pause (Time.Span.of_float 1.)
+          Time.pause (Time.Span.of_sec 1.)
         done
       end
   else
diff --git a/ocaml/Qmcchem_edit.ml b/ocaml/Qmcchem_edit.ml
index bcad9c7..8a757f3 100644
--- a/ocaml/Qmcchem_edit.ml
+++ b/ocaml/Qmcchem_edit.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 let file_header filename = Printf.sprintf
 "
diff --git a/ocaml/Qmcchem_forwarder.ml b/ocaml/Qmcchem_forwarder.ml
index 6a5719c..b5807d3 100644
--- a/ocaml/Qmcchem_forwarder.ml
+++ b/ocaml/Qmcchem_forwarder.ml
@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core
 
 let bind_socket ~socket_type ~socket ~address =
   let rec loop = function 
@@ -11,7 +11,7 @@ let bind_socket ~socket_type ~socket ~address =
         ZMQ.Socket.bind socket address; 
         loop (-1) 
       with 
-      | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_float 1. ; loop (i-1) ) 
+      | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_sec 1. ; loop (i-1) ) 
       | other_exception -> raise other_exception 
   in loop 10 
 
@@ -40,7 +40,7 @@ let run ezfio_filename dataserver =
   in
 
   (* Build qmc executable command *)
-  let prog, args = 
+  let prog, argv = 
       qmc, 
     [ qmc ; ezfio_filename ;
        Printf.sprintf "ipc://%s:%d" Qmcchem_config.dev_shm port ];
@@ -57,7 +57,7 @@ let run ezfio_filename dataserver =
     | Unix.Unix_error _ ->
         begin
           Unix.chdir tmpdir;
-          Time.pause @@ Time.Span.of_float 0.1;
+          Time.pause @@ Time.Span.of_sec 0.1;
           match (Sys.file_exists "PID") with
           | `No 
           | `Unknown -> ()
@@ -75,7 +75,7 @@ let run ezfio_filename dataserver =
                   begin
                     match Signal.send (Signal.of_system_int 0) (`Pid (Pid.of_int pid)) with
                     | `No_such_process ->  ()
-                    | _ -> ignore @@ Unix.exec ~prog ~args ()
+                    | _ -> ignore @@ Unix.exec ~prog ~argv ()
                   end
         end
   in
@@ -89,7 +89,7 @@ let run ezfio_filename dataserver =
 
   (* Fork a qmc *)
   ignore @@
-    Watchdog.fork_exec ~prog ~args ();
+    Watchdog.fork_exec ~prog ~argv ();
 
   (* If there are MICs, use them here (TODO) *)
   
diff --git a/ocaml/Qmcchem_info.ml b/ocaml/Qmcchem_info.ml
index 3a81458..a293791 100644
--- a/ocaml/Qmcchem_info.ml
+++ b/ocaml/Qmcchem_info.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 
 let run ezfio_filename = 
@@ -6,12 +6,12 @@ let run ezfio_filename =
   let qmcchem_info =
       Lazy.force Qmcchem_config.qmcchem_info
   in
-  let prog, args = 
+  let prog, argv = 
       qmcchem_info, 
     [ qmcchem_info ; ezfio_filename ]
   in
   ignore @@
-    Unix.exec ~prog ~args () 
+    Unix.exec ~prog ~argv () 
 
 
 let spec =
diff --git a/ocaml/Qmcchem_md5.ml b/ocaml/Qmcchem_md5.ml
index 54669a6..2f728cd 100644
--- a/ocaml/Qmcchem_md5.ml
+++ b/ocaml/Qmcchem_md5.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 let run ?c ?d ~l ~update ezfio_filename =
 
diff --git a/ocaml/Qmcchem_result.ml b/ocaml/Qmcchem_result.ml
index c1737d6..de93ec5 100644
--- a/ocaml/Qmcchem_result.ml
+++ b/ocaml/Qmcchem_result.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 (** Display a table that can be plotted by gnuplot *)
@@ -65,7 +65,7 @@ let display_cumulants ~range property =
   Printf.printf "Variance    = %16.10f\n" cum.(1);
   Printf.printf "Centered k3 = %16.10f\n" cum.(2);
   Printf.printf "Centered k4 = %16.10f\n" cum.(3);
-  print_newline ();
+  Printf.printf "\n%!";
   let n = 1. /. 12. *. cum.(2) *. cum.(2) +.
           1. /. 48. *. cum.(3) *. cum.(3)
   in
diff --git a/ocaml/Qmcchem_run.ml b/ocaml/Qmcchem_run.ml
index 03d9bc1..7a1a465 100644
--- a/ocaml/Qmcchem_run.ml
+++ b/ocaml/Qmcchem_run.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 let full_run ?(start_dataserver=true) ezfio_filename  = 
   (* Identify the job scheduler *)
@@ -36,13 +36,13 @@ let full_run ?(start_dataserver=true) ezfio_filename  =
 
 
       (* Start the data server *)
-      let prog, args = 
+      let prog, argv = 
         qmcchem,  [ qmcchem; "run" ; "-d" ; ezfio_filename]
       in
       let pid_dataserver = 
-        Watchdog.fork_exec ~prog ~args ()
+        Watchdog.fork_exec ~prog ~argv ()
       in
-      Printf.printf "%7d : %s\n%!" (Pid.to_int pid_dataserver) (String.concat ~sep:" " args)
+      Printf.printf "%7d : %s\n%!" (Pid.to_int pid_dataserver) (String.concat ~sep:" " argv)
     end;
   
 
@@ -83,7 +83,7 @@ let full_run ?(start_dataserver=true) ezfio_filename  =
   | n  -> 
     if (not (test_open_rep_socket ())) then
       begin
-       Time.pause (Time.Span.of_float 0.5);
+       Time.pause (Time.Span.of_sec 0.5);
        count (n-1);
       end
     else
@@ -94,7 +94,7 @@ let full_run ?(start_dataserver=true) ezfio_filename  =
      
 
   (* Start the qmc processes *)
-  let prog, args =
+  let prog, argv =
     let launcher =
       Launcher.(find () |> to_string)
     in
@@ -110,12 +110,12 @@ let full_run ?(start_dataserver=true) ezfio_filename  =
   in
   let pid_qmc = 
     try 
-      Watchdog.fork_exec ~prog ~args ()
+      Watchdog.fork_exec ~prog ~argv ()
     with
     | Unix.Unix_error _ ->
         begin
           let command = 
-            String.concat ~sep:" " args 
+            String.concat ~sep:" " argv 
           in
           Printf.printf "
 ============================================================
@@ -126,7 +126,7 @@ Error: Unable to run the following command
           Watchdog.kill () 
         end
   in
-  Printf.printf "%7d : %s\n%!" (Pid.to_int pid_qmc) (String.concat ~sep:" " args);  
+  Printf.printf "%7d : %s\n%!" (Pid.to_int pid_qmc) (String.concat ~sep:" " argv);  
 
   (* Wait for processes to finish *)
   Watchdog.join ()
diff --git a/ocaml/Qmcchem_stop.ml b/ocaml/Qmcchem_stop.ml
index 3224d29..2607eac 100644
--- a/ocaml/Qmcchem_stop.ml
+++ b/ocaml/Qmcchem_stop.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 
 let run ezfio_filename = 
diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml
index 19a5565..7ae4ffa 100644
--- a/ocaml/Qputils.ml
+++ b/ocaml/Qputils.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 let split_re = 
   Str.regexp " +"
diff --git a/ocaml/Random_variable.ml b/ocaml/Random_variable.ml
index c1ae8ad..b8e8a26 100644
--- a/ocaml/Random_variable.ml
+++ b/ocaml/Random_variable.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 open Qptypes
 
 type t = 
@@ -64,7 +64,7 @@ end = struct
       (x -. mu) *. ( x -. mu) /. sigma2
     in
     let pi = 
-      acos (-1.)
+      Float.acos (-1.)
     in
     let c = 
       1. /. (sqrt (sigma2 *. (pi +. pi)))
diff --git a/ocaml/Sample.ml b/ocaml/Sample.ml
index 14b709d..8f6818f 100644
--- a/ocaml/Sample.ml
+++ b/ocaml/Sample.ml
@@ -1,9 +1,9 @@
-open Core.Std
+open Core
 
 type t = 
   | One_dimensional  of float
   | Multidimensional of (float array * int)
-with sexp
+[@ deriving sexp]
 
 let dimension = function
   | One_dimensional _ -> 1
diff --git a/ocaml/Sample.mli b/ocaml/Sample.mli
index 27c965c..afa2308 100644
--- a/ocaml/Sample.mli
+++ b/ocaml/Sample.mli
@@ -1,4 +1,6 @@
-type t with sexp
+open Core
+
+type t [@@ deriving sexp]
 val to_float        : ?idx:int -> t -> float
 val to_float_array  : t -> float array
 val of_float        : float -> t
diff --git a/ocaml/Scheduler.ml b/ocaml/Scheduler.ml
index 06084c6..f20b184 100644
--- a/ocaml/Scheduler.ml
+++ b/ocaml/Scheduler.ml
@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core
 
 type t = 
   | SGE
diff --git a/ocaml/Watchdog.ml b/ocaml/Watchdog.ml
index 3a6a0aa..d0e92b8 100644
--- a/ocaml/Watchdog.ml
+++ b/ocaml/Watchdog.ml
@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core
 
 let _list    = ref [] ;;
 let _running = ref false;;
@@ -90,9 +90,9 @@ let del pid =
 ;;
 
 (** Fork and exec a new process *)
-let fork_exec ~prog ~args () =
+let fork_exec ~prog ~argv () =
   let pid = 
-    Unix.fork_exec ~prog ~args ()
+    Unix.fork_exec ~prog ~argv ()
   in
 
   let f () = 
diff --git a/ocaml/build.ninja b/ocaml/build.ninja
index 8238517..76786c0 100644
--- a/ocaml/build.ninja
+++ b/ocaml/build.ninja
@@ -1,7 +1,7 @@
 MAIN=qmcchem
 # Main program to build
 
-PACKAGES=-package core,cryptokit,str,ZMQ,sexplib.syntax
+PACKAGES=-package core,cryptokit,str,ZMQ
 #,ppx_sexp_conv
 # Required opam packages, for example:
 # PACKAGES=-package core,sexplib.syntax
@@ -10,7 +10,7 @@ THREAD=-thread
 # If you need threding support, use:
 # THREAD=-thread
 
-SYNTAX=-syntax camlp4o 
+SYNTAX=
 # If you need pre-processing, use:
 # SYNTAX=-syntax camlp4o
 
diff --git a/ocaml/ninja_ocaml.py b/ocaml/ninja_ocaml.py
index 00d992e..775b670 100755
--- a/ocaml/ninja_ocaml.py
+++ b/ocaml/ninja_ocaml.py
@@ -196,7 +196,7 @@ MAIN=
 
 PACKAGES=
 # Required opam packages, for example:
-# PACKAGES=-package core,sexplib.syntax
+# PACKAGES=-package core
 
 THREAD=
 # If you need threding support, use:
diff --git a/ocaml/qmcchem.ml b/ocaml/qmcchem.ml
index db1d35a..ad341b3 100644
--- a/ocaml/qmcchem.ml
+++ b/ocaml/qmcchem.ml
@@ -1,4 +1,4 @@
-open Core.Std
+open Core
 
 
 let command =
diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml
index cbcbd26..9889b92 100644
--- a/ocaml/qptypes_generator.ml
+++ b/ocaml/qptypes_generator.ml
@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core
 
 let input_data = "
 * Positive_float : float  
@@ -156,12 +156,12 @@ let untouched = "
 
 let template = format_of_string "
 module %s : sig
-  type t with sexp
+  type t [@@ deriving sexp]
   val to_%s : t -> %s
   val of_%s : %s %s -> t
   val to_string : t -> string
 end = struct
-  type t = %s with sexp
+  type t = %s [@@ deriving sexp]
   let to_%s x = x
   let of_%s %s x = ( %s x )
   let to_string x = %s.to_string x
@@ -199,13 +199,13 @@ let parse_input input=
 
 let ezfio_template = format_of_string "
 module %s : sig
-  type t with sexp
+  type t [@@ deriving sexp]
   val to_%s : t -> %s
   val get_max : unit -> %s
   val of_%s : ?min:%s -> ?max:%s -> %s -> t
   val to_string : t -> string
 end = struct
-  type t = %s with sexp
+  type t = %s [@@ deriving sexp]
   let to_string x = %s.to_string x
   let get_max () =
     if (Ezfio.has_%s ()) then
@@ -312,7 +312,7 @@ match msg with " ] @
 let () = 
   let input = 
     String.concat ~sep:"\n"
-    [ "open Core.Std\nlet warning = print_string\n\n" ;
+    [ "open Core\nlet warning = print_string\n\n" ;
       parse_input input_data ;
       parse_input_ezfio input_ezfio ;
       create_ezfio_handler ();
diff --git a/promela/qmcchem.pml b/promela/qmcchem.pml
new file mode 100644
index 0000000..a0c001b
--- /dev/null
+++ b/promela/qmcchem.pml
@@ -0,0 +1,271 @@
+#define NPROC 2
+#define BUFSIZE 4
+#define not_here false
+
+mtype = { NONE, TERMINATE, OK, TEST, ERROR, PROPERTY, WALKERS, EZFIO, GETWALKERS, REGISTER,
+          EZFIO_REPLY, UNREGISTER, STOPPING, STOPPED, QUEUED, RUNNING };
+
+typedef message_req {
+  mtype m = NONE;
+  byte value = 0;
+  chan reply = [BUFSIZE] of { mtype };
+}
+
+typedef message_pull {
+  mtype m = NONE;
+  byte value = 0;
+}
+
+chan dataserver_pull = [NPROC] of { message_pull };
+chan dataserver_req  = [NPROC] of { message_req };
+
+byte dataserver_status_pub;
+bit http_address = 0;
+bit killall_qmc = 0;
+bit killall_dataserver = 0;
+byte dataserver_status = QUEUED;
+byte dataserver_status_n_connected = 0;
+
+
+
+
+/* qmcchem process */
+active proctype qmcchem() {
+  byte reply = NONE;
+  byte dataserver_pid;
+  byte i,j;
+  message_req msg;
+
+  dataserver_pid = run dataserver();
+
+  /* Wait until ZMQ socket is open */
+  (http_address == 1);
+  do
+  :: (reply == OK) -> break
+  :: (reply == NONE) -> 
+        msg.m = TEST;
+        dataserver_req ! msg;
+        msg.reply ? reply ;
+        assert (reply == OK || reply == NONE)
+  od;
+  printf("Dataserver is ready.\n");
+
+  /* Start the QMC processes */
+
+  printf("qmcchem: Starting qmc processes.\n");
+  atomic {
+  i=0;
+  do
+  :: (i < NPROC) -> 
+    run qmc(); i++
+  :: else -> break
+  od;
+  }
+  printf("qmcchem: qmc processes started.\n");
+
+}
+
+
+
+
+
+
+
+/* dataserver process */
+proctype dataserver() {
+  byte reply = 0;
+  byte request = 0;
+  byte cont = 0;
+  byte reply_pid = 0;
+  message_req msg;
+
+  /* Simulate initialization */
+  http_address = 1;
+  dataserver_req ? msg;
+  msg.reply ! NONE ;
+
+  /* Status thread */
+  run dataserver_status_thread();
+  run dataserver_main_thread();
+}
+
+#define delay 5
+#define stop_time 100
+
+
+proctype dataserver_status_thread() {
+  byte count=0;
+  byte n_connected = 0;
+  byte time=0;
+
+  dataserver_status_pub = dataserver_status;
+  do
+  :: (dataserver_status == STOPPED) -> break
+  :: else  ->
+     time = (time < stop_time -> time+1 : time);
+     count++;
+     if
+     :: (count != delay) -> skip
+     :: else ->
+        count = 0;
+        if 
+        :: (dataserver_status == RUNNING &&
+            n_connected == dataserver_status_n_connected && 
+            time >= stop_time) ->
+            dataserver_status = STOPPING;
+            printf("Stop time reached : STOPPING\n")
+
+        :: (dataserver_status == STOPPING &&
+            n_connected != dataserver_status_n_connected && 
+            dataserver_status_n_connected == 0) ->
+            dataserver_status = STOPPED;
+            printf("No more connected clients : STOPPED\n")
+
+        :: (n_connected != dataserver_status_n_connected && 
+            dataserver_status_n_connected > 0) ->
+            n_connected = dataserver_status_n_connected;
+
+        :: else -> skip
+        fi
+     fi
+     dataserver_status_pub = dataserver_status;
+  od
+  printf ("End of dataserver_status_thread\n");
+}
+
+
+proctype dataserver_main_thread() {
+  byte time = 0;
+  mtype reply;
+  dataserver_status = QUEUED;
+  message_req msg;
+  message_pull pull;
+
+  /* Inform main process that the qmc processes can start (blocking recv) */
+  dataserver_req ? msg;
+  assert (msg.m == TEST);
+  msg.reply ! OK;
+
+  do
+  :: (dataserver_status == STOPPED && (!dataserver_pull ?[pull]) && (!dataserver_req ?[msg])) -> break
+  :: else ->
+    do
+    :: (dataserver_pull ?[pull]) -> 
+        dataserver_pull ? pull
+        printf("pull: "); printm(pull.m); printf("\n");
+        if
+        :: (pull.m == ERROR) -> skip;
+        :: (pull.m == WALKERS) -> skip
+        :: (pull.m == PROPERTY) -> skip;
+        fi
+    :: else -> break
+    od
+    if
+    :: (dataserver_req ?[msg]) -> 
+        dataserver_req ? msg;
+        printf("req : "); printm(msg.m); printf("\n");
+        if 
+        :: (msg.m == TEST) -> reply = OK
+        :: (msg.m == EZFIO) -> reply = EZFIO_REPLY
+        :: (msg.m == GETWALKERS) -> reply = WALKERS
+        :: (msg.m == REGISTER && dataserver_status == QUEUED ) ->
+              dataserver_status_n_connected++;
+              dataserver_status = RUNNING;
+              reply = OK;
+              printf("Status changed to RUNNING\n")
+        :: (msg.m == REGISTER && dataserver_status == RUNNING ) ->
+              dataserver_status_n_connected++;
+              reply = OK
+        :: (msg.m == REGISTER &&
+            (dataserver_status == STOPPED || dataserver_status == STOPPING) ) ->
+              dataserver_status_n_connected++; reply = ERROR;
+              printf("dataserver_req: register failed \n")
+        :: (msg.m == UNREGISTER) -> 
+            dataserver_status_n_connected--;
+            reply = OK;
+            if
+            :: (dataserver_status_n_connected == 0) ->
+                dataserver_status = STOPPED
+              printf("Status changed to STOPPED\n")
+            :: else -> skip
+            fi
+        :: else -> skip
+        fi;
+        msg.reply ! reply
+    :: else -> skip
+    fi
+  od
+}
+
+/* qmc processes */
+proctype qmc() {
+  byte status;
+  mtype reply;
+  message_req msg;
+  message_pull pull;
+
+  /* Init */
+  status = dataserver_status_pub;
+
+  msg.m = REGISTER;
+  dataserver_req ! msg;
+end:  msg.reply ? reply;
+  if 
+  :: (reply  == ERROR) -> goto exit;
+  :: else -> assert (reply == OK);
+  fi;
+
+  msg.m = EZFIO;
+  dataserver_req ! msg;
+  msg.reply ? reply;
+  if 
+  :: (reply  == ERROR) -> goto exit;
+  :: else -> assert (reply == EZFIO_REPLY);
+  fi;
+
+  msg.m = GETWALKERS;
+  dataserver_req ! msg;
+  msg.reply ? reply;
+  if 
+  :: (reply  == ERROR) -> goto exit;
+  :: else -> assert (reply == WALKERS);
+  fi;
+
+
+
+  /* Equilibration */
+  (dataserver_status_pub == RUNNING);
+
+  msg.m = EZFIO;
+  dataserver_req ! msg;
+  msg.reply ? reply;
+  if 
+  :: (reply  == ERROR) -> goto exit;
+  :: else -> assert (reply == EZFIO_REPLY);
+  fi;
+
+  status = dataserver_status_pub;
+
+  /* Cycles */
+  do
+  :: (status != RUNNING) -> break
+  :: else ->
+    pull.m = PROPERTY; pull.value = 0;
+    dataserver_pull ! pull;
+    pull.m = PROPERTY; pull.value =1 ;
+    dataserver_pull ! pull;
+    pull.m = WALKERS;
+    dataserver_pull ! pull;
+    status = dataserver_status_pub;
+  od;
+
+  /* Termination */
+  msg.m = UNREGISTER;
+  dataserver_req ! msg;
+  msg.reply ? reply;
+  assert (reply == OK);
+
+exit: skip
+}
+
+
diff --git a/scripts/create_properties_ezfio.py b/scripts/create_properties_ezfio.py
index 6872232..d2fda6c 100755
--- a/scripts/create_properties_ezfio.py
+++ b/scripts/create_properties_ezfio.py
@@ -67,7 +67,7 @@ file = open(tmp_filename,'w')
 # ----
 
 print >>file, """
-(* File generated by ${QMCCHEM_PATH}/src/create_properties.py. Do not
+(* File generated by ${QMCCHEM_PATH}/scripts/create_properties.py. Do not
    modify here
  *)
 
@@ -125,7 +125,7 @@ for p in properties:
 print >>file, """;;
 
 let of_string s = 
-  match (String.lowercase s) with 
+  match (String.lowercase_ascii s) with 
   | "cpu" -> Cpu
   | "wall" -> Wall
   | "accep" -> Accep"""
diff --git a/src/AO/ao_axis.irp.f b/src/AO/ao_axis.irp.f
index 9cfde0f..e46cbf2 100644
--- a/src/AO/ao_axis.irp.f
+++ b/src/AO/ao_axis.irp.f
@@ -45,7 +45,7 @@ subroutine pow_l(r,a,x1,x2,x3)
      x3 = 0.
      return
   end select
-end function
+end
 
 
  BEGIN_PROVIDER [ real, ao_axis_block, (ao_block_num_8) ]
diff --git a/src/JASTROW/jastrow_core.irp.f b/src/JASTROW/jastrow_core.irp.f
index f498a98..da8af7f 100644
--- a/src/JASTROW/jastrow_core.irp.f
+++ b/src/JASTROW/jastrow_core.irp.f
@@ -14,7 +14,7 @@
      jast_elec_Core_range(i) = 0.d0
    else
      double precision :: rc
-     double precision, parameter :: thresh = 0.5  ! function = thresh at rc
+     double precision, parameter :: thresh = 0.5d0  ! function = thresh at rc
      rc = min(0.8d0,max(4.0d0/nucl_charge(i), 0.25d0))
      jast_elec_Core_expo(i) = -1.d0/rc**2 * log(thresh)
      jast_elec_Core_range(i) = dsqrt(15.d0/jast_elec_Core_expo(i))
diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f
index 357b21d..12ebb73 100644
--- a/src/PROPERTIES/properties_energy.irp.f
+++ b/src/PROPERTIES/properties_energy.irp.f
@@ -270,21 +270,14 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ]
  BEGIN_DOC
  ! Zero-variance parameter on E_loc
  END_DOC
-   E_loc_zv  = E_loc + (E_trial-E_loc) * dmc_zv_weight
+   E_loc_zv = E_loc 
+   E_loc_zv += (E_trial-E_loc) * dmc_zv_weight
+!   E_loc_zv += - time_step*(E_trial**2 + 1.44341217940434 - E_loc**2)*dmc_zv_weight
 !  E_loc_zv(3)  = dmc_zv_weight_half
 !  E_loc_zv(:)  = 0.d0
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, E_loc_zv_diag ]
- implicit none
- BEGIN_DOC
- ! Zero-variance parameter on E_loc
- END_DOC
-  E_loc_zv_diag = E_trial
-
-END_PROVIDER
-
 
 
 
diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f
index 2764609..27397bc 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(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..c7e0e95 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(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
index 0efb370..484d7d8 100644
--- a/src/SAMPLING/pdmc_step.irp.f
+++ b/src/SAMPLING/pdmc_step.irp.f
@@ -109,9 +109,9 @@ END_SHELL
  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
+!  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)
 
@@ -234,13 +234,13 @@ 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)
+!    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
@@ -280,11 +280,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(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/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f
index 295d2ad..914cc8a 100644
--- a/src/SAMPLING/srmc_step.irp.f
+++ b/src/SAMPLING/srmc_step.irp.f
@@ -90,8 +90,8 @@ for p in properties:
  print  t.replace("$X",p[1])
 END_SHELL
 
- logical                        :: loop
- integer*8                      :: cpu0, cpu1, cpu2, count_rate, count_max
+ logical                      :: loop
+ integer*8                    :: cpu0, cpu1, cpu2, count_rate, count_max
 
  loop = .True.
  call system_clock(cpu0, count_rate, count_max)
@@ -320,11 +320,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(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
 
diff --git a/src/TOOLS/invert.irp.f b/src/TOOLS/invert.irp.f
index 98e1ff3..73d4b0a 100644
--- a/src/TOOLS/invert.irp.f
+++ b/src/TOOLS/invert.irp.f
@@ -1,3 +1,45 @@
+! This file contains the fast inversion routines of QMC=Chem for
+! small matrices. It may be downloaded here:
+! https://raw.githubusercontent.com/scemama/qmcchem/master/src/TOOLS/invert.irp.f
+! 
+! To use it in your Fortran code, you will need to~:
+! 1) rename it inverse.f90
+! 2) replace all $IRP_ALIGN occurences by 
+!    a) 16  for SSE4.2
+!    b) 32  for AVX or AVX2
+!    c) 64  for AVX-512
+!
+!
+! GPL license :
+! -------------
+!
+!   QMC=Chem : Quantum Monte Carlo for Chemistry
+!   Copyright (C) 2009 Anthony SCEMAMA 
+!
+!   This program is free software; you can redistribute it and/or modify
+!   it under the terms of the GNU General Public License as published by
+!   the Free Software Foundation; either version 2 of the License, or
+!   (at your option) any later version.
+!
+!   This program is distributed in the hope that it will be useful,
+!   but WITHOUT ANY WARRANTY; without even the implied warranty of
+!   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+!   GNU General Public License for more details.
+!
+!   You should have received a copy of the GNU General Public License along
+!   with this program; if not, write to the Free Software Foundation, Inc.,
+!   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!   Anthony Scemama
+!   LCPQ - IRSAMC - CNRS
+!   Universite Paul Sabatier
+!   118, route de Narbonne      
+!   31062 Toulouse Cedex 4      
+!   scemama@irsamc.ups-tlse.fr
+!
+!
+
+
 subroutine invert(a,LDA,na,det_l)
  implicit none
  double precision, intent(inout) :: a (LDA,na)
@@ -43,7 +85,6 @@ subroutine invert_general(a,LDA,na,det_l)
  integer          :: ipiv(LDA)
 !DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: ipiv
   integer          :: lwork
-  double precision :: f
   integer          :: i, j
   call dgetrf(na, na, a, LDA, ipiv, inf )
   det_l = 1.d0
@@ -74,7 +115,6 @@ subroutine sinvert(a,LDA,na,det_l)
  integer          :: ipiv(LDA)
 !DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: ipiv
   integer          :: lwork
-  real             :: f
   integer          :: i, j
   call sgetrf(na, na, a, LDA, ipiv, inf )
   det_l = 1.d0
@@ -113,7 +153,6 @@ subroutine invert2(a,LDA,na,det_l)
  double precision :: det_l
  double precision :: b(2,2)
 
- double precision :: f
  b(1,1) = a(1,1)
  b(2,1) = a(2,1)
  b(1,2) = a(1,2)
@@ -134,7 +173,6 @@ subroutine invert3(a,LDA,na,det_l)
  double precision :: b(4,3)
  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b
  integer :: i
- double precision :: f
  det_l = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) &
         -a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) &
         +a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
@@ -166,7 +204,6 @@ subroutine invert4(a,LDA,na,det_l)
  double precision :: b(4,4)
  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b
  integer :: i,j
- double precision :: f
  det_l =  a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))  &
                  -a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))  &
                  +a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))) &
@@ -217,7 +254,6 @@ subroutine invert5(a,LDA,na,det_l)
  double precision :: b(5,5)
  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b
  integer :: i,j
- double precision :: f
  det_l = a(1,1)*(a(2,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( &
  a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))- &
  a(2,3)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)- &
@@ -410,7 +446,6 @@ subroutine invert_update(a,LDA,na,det_l,b)
  integer          :: ipiv(LDA)
 !DIR$ ATTRIBUTES ALIGN: $IRP_ALIGN :: ipiv
   integer          :: lwork
-  double precision :: f
   integer          :: i, j
   double precision :: bold(LDA,na)
   double precision :: ba(LDA,na)
diff --git a/src/det.irp.f b/src/det.irp.f
index 0982d07..9b1f6d0 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -1114,6 +1114,9 @@ end
  endif
  !DIR$ FORCEINLINE
  call bitstring_to_list ( psi_det_alpha(1,det_i), mo_list_alpha_curr, l, N_int )
+ if (l /= elec_alpha_num) then
+   stop 'error in number of alpha electrons'
+ endif
 
 END_PROVIDER
 
@@ -1132,8 +1135,12 @@ END_PROVIDER
  else
    mo_list_beta_prev = 0
  endif
+
  !DIR$ FORCEINLINE
  call bitstring_to_list ( psi_det_beta(1,det_j), mo_list_beta_curr, l, N_int )
+ if (l /= elec_beta_num) then
+   stop 'error in number of beta electrons'
+ endif
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, det_alpha_value_curr ]
@@ -1401,15 +1408,10 @@ END_PROVIDER
     endif
     
     det_alpha_value(det_i) = det_alpha_value_curr
-    do i=1,elec_alpha_num
-      !DIR$ VECTOR ALIGNED
-      do k=1,4
-        det_alpha_grad_lapl(k,i,det_i) = det_alpha_grad_lapl_curr(k,i)
-      enddo
-      if (do_pseudo) then
-        det_alpha_pseudo(i,det_i) = det_alpha_pseudo_curr(i) 
-      endif
-    enddo
+    det_alpha_grad_lapl(:,:,det_i) = det_alpha_grad_lapl_curr(:,:)
+    if (do_pseudo) then
+      det_alpha_pseudo(:,det_i) = det_alpha_pseudo_curr(:) 
+    endif
    
   enddo
   
@@ -1453,15 +1455,10 @@ END_PROVIDER
     endif
     
     det_beta_value(det_j)  = det_beta_value_curr
-    !DIR$ LOOP COUNT (200)
-    do i=elec_alpha_num+1,elec_num
-      do k=1,4
-        det_beta_grad_lapl(k,i,det_j) = det_beta_grad_lapl_curr(k,i)
-      enddo
-      if (do_pseudo) then
-        det_beta_pseudo(i,det_j) = det_beta_pseudo_curr(i) 
-      endif
-    enddo
+    det_beta_grad_lapl(:,:,det_j) = det_beta_grad_lapl_curr(:,:)
+    if (do_pseudo) then
+      det_beta_pseudo(:,det_j) = det_beta_pseudo_curr(:) 
+    endif
 
   enddo
   
@@ -1536,53 +1533,64 @@ END_PROVIDER
   DaC = 0.d0
   CDb = 0.d0
 
-  det_num4 = iand(det_num,not(3))
-  !DIR$ VECTOR ALIGNED
-  do k=1,det_num4,4
-    i1 = det_coef_matrix_rows(k  )
-    i2 = det_coef_matrix_rows(k+1)
-    i3 = det_coef_matrix_rows(k+2)
-    i4 = det_coef_matrix_rows(k+3)
-    j1 = det_coef_matrix_columns(k  )
-    j2 = det_coef_matrix_columns(k+1)
-    j3 = det_coef_matrix_columns(k+2)
-    j4 = det_coef_matrix_columns(k+3)
-    if ( (j1 == j2).and.(j1 == j3).and.(j1 == j4) ) then
-      f = det_beta_value (j1)
-      CDb(i1) = CDb(i1) + det_coef_matrix_values(k  )*f
-      CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*f
-      CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*f
-      CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*f
+  if (det_num < ishft(det_alpha_num*det_beta_num,-2)) then
 
-      if ( ((i2-i1) == 1).and.((i3-i1) == 2).and.((i4-i1) == 3) ) then
-        DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
-         + det_coef_matrix_values(k+1)*det_alpha_value(i1+1) &
-         + det_coef_matrix_values(k+2)*det_alpha_value(i1+2) &
-         + det_coef_matrix_values(k+3)*det_alpha_value(i1+3)
+    det_num4 = iand(det_num,not(3))
+    !DIR$ VECTOR ALIGNED
+    do k=1,det_num4,4
+      i1 = det_coef_matrix_rows(k  )
+      i2 = det_coef_matrix_rows(k+1)
+      i3 = det_coef_matrix_rows(k+2)
+      i4 = det_coef_matrix_rows(k+3)
+      j1 = det_coef_matrix_columns(k  )
+      j2 = det_coef_matrix_columns(k+1)
+      j3 = det_coef_matrix_columns(k+2)
+      j4 = det_coef_matrix_columns(k+3)
+      if ( (j1 == j2).and.(j1 == j3).and.(j1 == j4) ) then
+        f = det_beta_value (j1)
+        CDb(i1) = CDb(i1) + det_coef_matrix_values(k  )*f
+        CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*f
+        CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*f
+        CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*f
+
+        if ( ((i2-i1) == 1).and.((i3-i1) == 2).and.((i4-i1) == 3) ) then
+          DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
+          + det_coef_matrix_values(k+1)*det_alpha_value(i1+1) &
+          + det_coef_matrix_values(k+2)*det_alpha_value(i1+2) &
+          + det_coef_matrix_values(k+3)*det_alpha_value(i1+3)
+        else
+          DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
+          + det_coef_matrix_values(k+1)*det_alpha_value(i2) &
+          + det_coef_matrix_values(k+2)*det_alpha_value(i3) &
+          + det_coef_matrix_values(k+3)*det_alpha_value(i4)
+        endif
       else
-        DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
-         + det_coef_matrix_values(k+1)*det_alpha_value(i2) &
-         + det_coef_matrix_values(k+2)*det_alpha_value(i3) &
-         + det_coef_matrix_values(k+3)*det_alpha_value(i4)
+        DaC(j1) = DaC(j1) + det_coef_matrix_values(k  )*det_alpha_value(i1)
+        DaC(j2) = DaC(j2) + det_coef_matrix_values(k+1)*det_alpha_value(i2)
+        DaC(j3) = DaC(j3) + det_coef_matrix_values(k+2)*det_alpha_value(i3)
+        DaC(j4) = DaC(j4) + det_coef_matrix_values(k+3)*det_alpha_value(i4)
+        CDb(i1) = CDb(i1) + det_coef_matrix_values(k  )*det_beta_value (j1)
+        CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*det_beta_value (j2)
+        CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*det_beta_value (j3)
+        CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*det_beta_value (j4)
       endif
-    else
-      DaC(j1) = DaC(j1) + det_coef_matrix_values(k  )*det_alpha_value(i1)
-      DaC(j2) = DaC(j2) + det_coef_matrix_values(k+1)*det_alpha_value(i2)
-      DaC(j3) = DaC(j3) + det_coef_matrix_values(k+2)*det_alpha_value(i3)
-      DaC(j4) = DaC(j4) + det_coef_matrix_values(k+3)*det_alpha_value(i4)
-      CDb(i1) = CDb(i1) + det_coef_matrix_values(k  )*det_beta_value (j1)
-      CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*det_beta_value (j2)
-      CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*det_beta_value (j3)
-      CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*det_beta_value (j4)
-    endif
-  enddo
+    enddo
 
-  do k=det_num4+1,det_num
-    i = det_coef_matrix_rows(k)
-    j = det_coef_matrix_columns(k)
-    DaC(j) = DaC(j) + det_coef_matrix_values(k)*det_alpha_value(i)
-    CDb(i) = CDb(i) + det_coef_matrix_values(k)*det_beta_value (j)
-  enddo
+    do k=det_num4+1,det_num
+      i = det_coef_matrix_rows(k)
+      j = det_coef_matrix_columns(k)
+      DaC(j) = DaC(j) + det_coef_matrix_values(k)*det_alpha_value(i)
+      CDb(i) = CDb(i) + det_coef_matrix_values(k)*det_beta_value (j)
+    enddo
+
+  else
+
+    call dgemv('T',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, &
+      size(det_coef_matrix_dense,1), det_alpha_value, 1, 0.d0, DaC, 1)
+    call dgemv('N',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, &
+      size(det_coef_matrix_dense,1), det_beta_value, 1, 0.d0, CDb, 1)
+
+  endif
 
   ! Value
   ! -----
diff --git a/src/mo.irp.f b/src/mo.irp.f
index eb300af..4b9bba7 100644
--- a/src/mo.irp.f
+++ b/src/mo.irp.f
@@ -701,13 +701,13 @@ subroutine sparse_full_mv(A,LDA,                                     &
   
   ! LDC and LDA have to be factors of simd_sp
   
-  IRP_IF NO_PREFETCH
-  IRP_ELSE
-    call MM_PREFETCH (A(j,indices(1)),3)
-    call MM_PREFETCH (A(j,indices(2)),3)
-    call MM_PREFETCH (A(j,indices(3)),3)
-    call MM_PREFETCH (A(j,indices(4)),3)
-  IRP_ENDIF
+!  IRP_IF NO_PREFETCH
+!  IRP_ELSE
+!    call MM_PREFETCH (A(1,indices(1)),3)
+!    call MM_PREFETCH (A(1,indices(2)),3)
+!    call MM_PREFETCH (A(1,indices(3)),3)
+!    call MM_PREFETCH (A(1,indices(4)),3)
+!  IRP_ENDIF
 
   !DIR$ SIMD
   do j=1,LDC
@@ -757,13 +757,13 @@ subroutine sparse_full_mv(A,LDA,                                     &
       !DIR$ VECTOR ALIGNED
       !DIR$ SIMD FIRSTPRIVATE(d11,d21,d31,d41)
       do j=1,$IRP_ALIGN/4
-        IRP_IF NO_PREFETCH
-        IRP_ELSE
-          call MM_PREFETCH (A(j+k,indices(kao+4)),3)
-          call MM_PREFETCH (A(j+k,indices(kao+5)),3)
-          call MM_PREFETCH (A(j+k,indices(kao+6)),3)
-          call MM_PREFETCH (A(j+k,indices(kao+7)),3)
-        IRP_ENDIF
+!        IRP_IF NO_PREFETCH
+!        IRP_ELSE
+!          call MM_PREFETCH (A(j+k,indices(kao+4)),3)
+!          call MM_PREFETCH (A(j+k,indices(kao+5)),3)
+!          call MM_PREFETCH (A(j+k,indices(kao+6)),3)
+!          call MM_PREFETCH (A(j+k,indices(kao+7)),3)
+!        IRP_ENDIF
         C1(j+k) = C1(j+k) + A(j+k,k_vec(1))*d11 + A(j+k,k_vec(2))*d21&
             +   A(j+k,k_vec(3))*d31 + A(j+k,k_vec(4))*d41
       enddo
diff --git a/src/simulation.irp.f b/src/simulation.irp.f
index d3ba571..88b74d0 100644
--- a/src/simulation.irp.f
+++ b/src/simulation.irp.f
@@ -264,7 +264,7 @@ END_PROVIDER
  nucl_fitcusp_factor = 0.
  call get_simulation_nucl_fitcusp_factor(nucl_fitcusp_factor)
  do_nucl_fitcusp = nucl_fitcusp_factor > 0.
- call info(irp_here,'nucl_fitcusp_factor',nucl_fitcusp_factor)
+ call rinfo(irp_here,'nucl_fitcusp_factor',nucl_fitcusp_factor)
 
 END_PROVIDER
 
diff --git a/src/wf.irp.f b/src/wf.irp.f
index cee557b..5e117f5 100644
--- a/src/wf.irp.f
+++ b/src/wf.irp.f
@@ -80,6 +80,21 @@ END_PROVIDER
   deallocate(buffer)
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, det_coef_matrix_dense, (det_alpha_num, det_beta_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Dense version of det_coef_matrix
+ END_DOC
+ integer :: i,j,k
+ det_coef_matrix_dense = 0.d0
+ do k=1,det_num
+   i = det_coef_matrix_rows(k)
+   j = det_coef_matrix_columns(k)
+   det_coef_matrix_dense(i,j) = det_coef_matrix_values(k)
+ enddo
+END_PROVIDER
+
+
 BEGIN_PROVIDER [ integer, det_num ]
   implicit none
   BEGIN_DOC