From 55bf56760ac2356ae92446658618afa0a2e78896 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Sat, 18 Jun 2016 01:05:55 +0200
Subject: [PATCH 01/17] Bug in gcc version with GCC 5 corrected

---
 install/scripts/install_ocaml.sh | 40 ++++++++++++++++++++++++++++++--
 1 file changed, 38 insertions(+), 2 deletions(-)

diff --git a/install/scripts/install_ocaml.sh b/install/scripts/install_ocaml.sh
index 5718ed3..78e1ef9 100755
--- a/install/scripts/install_ocaml.sh
+++ b/install/scripts/install_ocaml.sh
@@ -7,8 +7,44 @@ 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 2
+        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

From beae2ab7bf2cad8c1d7fe37a1b9b59324708739d Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Sat, 18 Jun 2016 22:40:22 +0200
Subject: [PATCH 02/17] GCC 5.3.1

---
 install/scripts/install_ocaml.sh | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/install/scripts/install_ocaml.sh b/install/scripts/install_ocaml.sh
index 78e1ef9..e6f0f0b 100755
--- a/install/scripts/install_ocaml.sh
+++ b/install/scripts/install_ocaml.sh
@@ -6,8 +6,6 @@ set -e
 cd .. ; QMCCHEM_PATH="$PWD" ; cd -
 PACKAGES="core cryptokit ocamlfind sexplib" # ppx_sexp_conv"
 
-declare -i i
-
 # return 0 if program version is equal or greater than check version
 check_version () {
     if [[ $1 == $2 ]]
@@ -34,7 +32,7 @@ check_version () {
         fi
         if ((10#${ver1[i]} < 10#${ver2[i]}))
         then
-            return 2
+            return 0
         fi
     done
     return 0

From b577984de6a0b5a962100269591d9df636486b5b Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 3 Oct 2016 00:34:42 +0200
Subject: [PATCH 03/17] Updated ZMQ

---
 install/build.ninja            |  7 +++----
 install/scripts/install_zmq.sh | 10 +++++-----
 2 files changed, 8 insertions(+), 9 deletions(-)

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_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

From dcf19db6d625743dfc7316b537b4c1ddc880a634 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 28 Nov 2016 15:08:55 +0100
Subject: [PATCH 04/17] Commented prefetch instructions

---
 src/det.irp.f |  7 +++++++
 src/mo.irp.f  | 28 ++++++++++++++--------------
 2 files changed, 21 insertions(+), 14 deletions(-)

diff --git a/src/det.irp.f b/src/det.irp.f
index 0982d07..85494f2 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 ]
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

From 4c42654401513d709721f90ea37c144f52007e9b Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 28 Nov 2016 19:45:09 +0100
Subject: [PATCH 05/17] Cleaned PDMC

---
 src/PROPERTIES/properties_energy.irp.f |  4 +++-
 src/SAMPLING/pdmc_step.irp.f           | 20 ++++++++++----------
 2 files changed, 13 insertions(+), 11 deletions(-)

diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f
index 357b21d..7f42053 100644
--- a/src/PROPERTIES/properties_energy.irp.f
+++ b/src/PROPERTIES/properties_energy.irp.f
@@ -270,7 +270,9 @@ 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
 
diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f
index 0efb370..9c9b4a8 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

From bb4b75e86b2fd1130d8dc069b02a8a84d24a7983 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 28 Dec 2016 16:57:53 +0100
Subject: [PATCH 06/17] Fixed Bug when block is too long

---
 src/SAMPLING/dmc_step.irp.f  |   5 +-
 src/SAMPLING/fkmc_step.irp.f |   4 +-
 src/SAMPLING/pdmc_step.irp.f | 380 +++++++++++++++++++++++++++++++++++
 src/SAMPLING/srmc_step.irp.f |   4 +-
 src/SAMPLING/vmc_step.irp.f  |   4 +-
 5 files changed, 389 insertions(+), 8 deletions(-)
 create mode 100644 src/SAMPLING/pdmc_step.irp.f

diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f
index 2764609..7044080 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(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
 
+
   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..3f57d20 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(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/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f
new file mode 100644
index 0000000..e4c3e70
--- /dev/null
+++ b/src/SAMPLING/pdmc_step.irp.f
@@ -0,0 +1,380 @@
+! Providers of *_pdmc_block_walk
+!==============================
+BEGIN_SHELL [ /usr/bin/python ]
+from properties import *
+
+t = """
+ BEGIN_PROVIDER [ $T, $X_pdmc_block_walk $D1 ]
+&BEGIN_PROVIDER [ $T, $X_pdmc_block_walk_kahan $D2 ]
+&BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk $D1 ]
+&BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk_kahan $D2 ]
+ implicit none
+ BEGIN_DOC  
+! pdMC averages of $X. Computed in E_loc_pdmc_block_walk
+ END_DOC
+ $X_pdmc_block_walk = 0.d0
+ $X_pdmc_block_walk_kahan = 0.d0
+ $X_2_pdmc_block_walk = 0.d0
+ $X_2_pdmc_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_pdmc_block_walk       ]
+&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk     ]
+&BEGIN_PROVIDER [ double precision, E_loc_pdmc_block_walk_kahan, (3) ]
+&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk_kahan, (3) ]
+ implicit none
+ include '../types.F'
+ BEGIN_DOC
+! Properties averaged over the block using the PDMC method
+ END_DOC
+
+  real, allocatable :: elec_coord_tmp(:,:,:)
+  integer :: mod_align
+  double precision :: E_loc_save(4,walk_num_dmc_max)
+  double precision :: psi_value_save(walk_num)
+  double precision :: psi_value_save_tmp(walk_num)
+  double precision :: pdmc_weight(walk_num)
+  double precision, allocatable :: psi_grad_psi_inv_save(:,:,:)
+  double precision, allocatable :: psi_grad_psi_inv_save_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 :: psi_value_save
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save_tmp
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: pdmc_weight
+  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) )
+  psi_value_save = 0.d0
+  psi_value_save_tmp = 0.d0
+  pdmc_weight = 1.d0
+
+! Initialization
+ if (vmc_algo /= t_Brownian) then
+   call abrt(irp_here,'PDMC 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_pdmc_block_walk = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_pdmc_block_walk_kahan = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_2_pdmc_block_walk = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_2_pdmc_block_walk_kahan = 0.d0
+ endif
+"""
+for p in properties:
+ print  t.replace("$X",p[1])
+END_SHELL
+
+ logical                        :: loop
+ integer*8                      :: cpu0, cpu1, cpu2, count_rate, count_max
+
+ loop = .True.
+ call system_clock(cpu0, count_rate, count_max)
+ cpu2 = cpu0
+
+ block_weight = 0.d0
+
+ real, external                 :: accep_rate
+ double precision               :: delta, thr
+
+ thr = 2.d0/time_step_sq
+
+ logical :: first_loop
+ first_loop = .True.
+ if (walk_num > 1) then
+   call abrt(irp_here,'walk_num > 1')
+ 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
+
+ do while (loop)
+
+  i_walk = 1
+    
+    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(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
+      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
+      psi_value_save(i_walk) = psi_value
+      E_loc_save(:,i_walk) = E_loc
+    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
+
+!2    delta = (E_loc+E_loc_save(1,i_walk))*0.5d0
+!3    delta = (5.d0 * E_loc + 8.d0 * E_loc_save(1,i_walk) - E_loc_save(2,i_walk))/12.d0
+      delta = (9.d0*E_loc+19.d0*E_loc_save(1,i_walk)-5.d0*E_loc_save(2,i_walk)+E_loc_save(3,i_walk))/24.d0
+!      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
+       pdmc_weight(i_walk) = dexp(-dtime_step*delta)
+     else
+       pdmc_weight(i_walk) = 2.d0-dexp(dtime_step*delta)
+     endif
+     elec_coord(elec_num+1,1) += p*time_step
+     elec_coord(elec_num+1,2)  = E_loc
+     elec_coord(elec_num+1,3)  = pdmc_weight(i_walk) * pdmc_pop_weight_mult(pdmc_n_diag)
+     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(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
+
+     if (dabs(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag)) > 1.d-15) then
+       dmc_zv_weight = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag))
+       dmc_zv_weight_half = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag/2))
+     else
+       dmc_zv_weight = 0.d0
+       dmc_zv_weight_half = 0.d0
+     endif
+     TOUCH dmc_zv_weight dmc_zv_weight_half
+
+!    do i=1,pdmc_n_diag+1
+!      E_loc_zv(i) =   E_loc * pdmc_pop_weight_mult(i-1) * pdmc_weight(i_walk) * dmc_zv_weight + (E_trial-E_loc) * dmc_zv_weight
+!      E_loc_zv(i+pdmc_n_diag+1) = pdmc_pop_weight_mult(i-1) * pdmc_weight(i_walk) * dmc_zv_weight
+!    enddo
+
+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_pdmc_block_walk    += $X * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk)
+   !  $X_2_pdmc_block_walk  += $X_2 * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk)
+   ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
+   
+      $X_pdmc_block_walk_kahan($D2 3) = $X * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) - $X_pdmc_block_walk_kahan($D2 1)
+      $X_pdmc_block_walk_kahan($D2 2) = $X_pdmc_block_walk $D1  + $X_pdmc_block_walk_kahan($D2 3)
+      $X_pdmc_block_walk_kahan($D2 1) = ($X_pdmc_block_walk_kahan($D2 2) - $X_pdmc_block_walk $D1 ) &
+          - $X_pdmc_block_walk_kahan($D2 3)
+      $X_pdmc_block_walk $D1  =  $X_pdmc_block_walk_kahan($D2 2) 
+   
+   
+      $X_2_pdmc_block_walk_kahan($D2 3) = $X_2 * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) - $X_2_pdmc_block_walk_kahan($D2 1)
+      $X_2_pdmc_block_walk_kahan($D2 2) = $X_2_pdmc_block_walk $D1 + $X_2_pdmc_block_walk_kahan($D2 3)
+      $X_2_pdmc_block_walk_kahan($D2 1) = ($X_2_pdmc_block_walk_kahan($D2 2) - $X_2_pdmc_block_walk $D1 ) &
+          - $X_2_pdmc_block_walk_kahan($D2 3)
+      $X_2_pdmc_block_walk $D1 =  $X_2_pdmc_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 += 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)
+
+!  else
+!     pdmc_weight(i_walk) = 1.d0
+!     pdmc_pop_weight(:,:) = 1.d0
+!     pdmc_pop_weight_mult(:) = 1.d0
+!  endif
+   
+
+  do k=1,pdmc_n_diag
+    ! Move to the next projection step
+    if (pdmc_projection(pdmc_n_diag) > 0) then
+      pdmc_projection_step(k) = mod(pdmc_projection_step(k),pdmc_projection(k))+1
+    else
+      pdmc_projection_step(k) = 1
+    endif
+
+    ! Eventually, recompute the weight of the population
+    if (pdmc_projection_step(k) == k) then
+      pdmc_pop_weight_mult(k) = 1.d0
+      do l=1,pdmc_projection(k)
+        pdmc_pop_weight_mult(k) *= pdmc_pop_weight(l,k)
+      enddo
+    endif
+    ! Remove contribution of the old value of the weight at the new
+    ! projection step
+    
+    pdmc_pop_weight_mult(k) *= 1.d0/pdmc_pop_weight(pdmc_projection_step(k),k)
+    pdmc_pop_weight(pdmc_projection_step(k),k) = pdmc_weight(i_walk)/dble(walk_num)
+
+    ! Update the running population weight
+    pdmc_pop_weight_mult(k) *= pdmc_pop_weight(pdmc_projection_step(k),k)
+
+  enddo
+
+
+  call system_clock(cpu1, count_rate, count_max)
+  if (cpu1 < cpu0) then
+    cpu1 = cpu1+cpu0
+  endif
+  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 = loop.and.(do_run == t_Running)
+    cpu2 = cpu1
+  endif
+
+  SOFT_TOUCH elec_coord_full pdmc_pop_weight_mult 
+
+  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_pdmc_block_walk   *= factor
+   $X_2_pdmc_block_walk *= factor
+ endif
+"""
+for p in properties:
+ print  t.replace("$X",p[1])
+END_SHELL
+
+!  H(0,0) = H(3,3)
+!  H(1,0) = H(4,3)
+!  H(0,1) = H(3,4)
+!  H(1,1) = H(4,4)
+!  S(0,0) = S(3,3)
+!  S(1,0) = S(4,3)
+!  S(0,1) = S(3,4)
+!  S(1,1) = S(4,4)
+!
+!  print *,  H(0,0)/S(0,0)
+!  print *,  H(1,1)/S(1,1)
+!  print *,  ''
+!
+!  call dsygv(1, 'N', 'U', pdmc_n_diag/2+1, H, pdmc_n_diag/2+1, S, pdmc_n_diag/2+1, w, work, 3*(pdmc_n_diag+1), info)
+!  call dsygv(1, 'N', 'U', 2, H, pdmc_n_diag/2+1, S, pdmc_n_diag/2+1, w, work, 3*(pdmc_n_diag+1), info)
+!  E_loc_zv_diag_pdmc_block_walk = w(0)
+!  print *,  w
+
+ deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp )
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer, pdmc_projection, (pdmc_n_diag) ]
+&BEGIN_PROVIDER [ integer, pdmc_projection_step, (pdmc_n_diag) ]
+ implicit none
+ BEGIN_DOC  
+! Number of projection steps for PDMC
+ END_DOC
+ real :: pdmc_projection_time
+ pdmc_projection_time = 1.
+ call get_simulation_srmc_projection_time(pdmc_projection_time)
+ pdmc_projection(pdmc_n_diag) = int( pdmc_projection_time/time_step)
+ integer :: k
+ do k=1,pdmc_n_diag-1
+   pdmc_projection(k) = k*pdmc_projection(pdmc_n_diag)/pdmc_n_diag
+ enddo
+ pdmc_projection_step(:) = 0
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection(pdmc_n_diag)+1,pdmc_n_diag) ]
+ implicit none
+ BEGIN_DOC  
+! Population weight of PDMC
+ END_DOC
+ pdmc_pop_weight(:,:) = 1.d0
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, pdmc_pop_weight_mult, (0:pdmc_n_diag) ]
+ implicit none
+ BEGIN_DOC  
+! Population weight of PDMC
+ END_DOC
+ pdmc_pop_weight_mult(:) = 1.d0
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, pdmc_n_diag ]
+ implicit none
+ BEGIN_DOC
+! Size of the matrix to diagonalize
+ END_DOC
+ pdmc_n_diag = 8
+END_PROVIDER
+
+
+
diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f
index 836ddae..1300fe2 100644
--- a/src/SAMPLING/srmc_step.irp.f
+++ b/src/SAMPLING/srmc_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(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/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
 

From a1192beb8b8672bd84c17ed4d376484eb0194cda Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 29 Dec 2016 01:43:34 +0100
Subject: [PATCH 07/17] Promela model

---
 promela/qmcchem.pml | 271 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 271 insertions(+)
 create mode 100644 promela/qmcchem.pml

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
+}
+
+

From c59606bb36863dfa59309081f6b71bac943f95b6 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 29 Dec 2016 01:48:59 +0100
Subject: [PATCH 08/17] Removed useless property

---
 src/PROPERTIES/properties_energy.irp.f | 9 ---------
 1 file changed, 9 deletions(-)

diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f
index 7f42053..12ebb73 100644
--- a/src/PROPERTIES/properties_energy.irp.f
+++ b/src/PROPERTIES/properties_energy.irp.f
@@ -278,15 +278,6 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ]
 
 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
-
 
 
 

From d783b19cc8b1a14322e9a195e2f530d774740a79 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Sun, 26 Mar 2017 21:36:23 +0200
Subject: [PATCH 09/17] ZMQ versio update

---
 install/build.ninja                  |  4 ++--
 install/scripts/install_ocaml_zmq.sh |  2 +-
 install/scripts/install_zmq.sh       | 11 ++++++-----
 src/JASTROW/jastrow_core.irp.f       |  2 +-
 4 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/install/build.ninja b/install/build.ninja
index 351db7d..758982d 100644
--- a/install/build.ninja
+++ b/install/build.ninja
@@ -7,8 +7,8 @@ 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_ZMQ   ="http://download.zeromq.org/zeromq-4.0.7.tar.gz"
+URL_ZMQ   ="http://download.zeromq.org/zeromq-4.1.4.tar.gz"
 URL_F77ZMQ="https://github.com/scemama/f77_zmq/archive/v4.1.3.tar.gz"
 
 # Rules
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..c3b99ea 100755
--- a/install/scripts/install_zmq.sh
+++ b/install/scripts/install_zmq.sh
@@ -3,7 +3,8 @@
 TARGET=zmq
 function _install()
 {
-  LIBVERSION=4
+#  LIBVERSION=4
+  LIBVERSION=5
   cd .. ; QMCCHEM_PATH="$PWD" ; cd -
   set +u
   export C_INCLUDE_PATH="${C_INCLUDE_PATH}":./
@@ -14,10 +15,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/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))

From fed153e513f0e96b33f7e51049e87378e75eeb9f Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 14 Jun 2017 01:11:26 +0200
Subject: [PATCH 10/17] New IRPF90

---
 src/AO/ao_axis.irp.f | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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) ]

From a89daa602be30bf127b00f56fcfb54988b530658 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 18 Jul 2017 18:08:34 +0200
Subject: [PATCH 11/17] Added header to inverse

---
 src/TOOLS/invert.irp.f | 42 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 42 insertions(+)

diff --git a/src/TOOLS/invert.irp.f b/src/TOOLS/invert.irp.f
index 98e1ff3..202c0cb 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)

From 0e2ee815f69e3732361ec0c4886a7392cca10efa Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 5 Sep 2017 17:39:11 +0200
Subject: [PATCH 12/17] Unused variable

---
 src/TOOLS/invert.irp.f | 7 -------
 1 file changed, 7 deletions(-)

diff --git a/src/TOOLS/invert.irp.f b/src/TOOLS/invert.irp.f
index 98e1ff3..0c8e3a3 100644
--- a/src/TOOLS/invert.irp.f
+++ b/src/TOOLS/invert.irp.f
@@ -43,7 +43,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 +73,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 +111,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 +131,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 +162,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 +212,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 +404,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)

From 9d186b37596531f4bf52a4e14aab0907de8b475d Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 9 Oct 2017 21:39:00 +0200
Subject: [PATCH 13/17] DGEMV for dense matrices

---
 src/det.irp.f | 97 ++++++++++++++++++++++++++++-----------------------
 src/wf.irp.f  | 15 ++++++++
 2 files changed, 69 insertions(+), 43 deletions(-)

diff --git a/src/det.irp.f b/src/det.irp.f
index 85494f2..17410f0 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -1543,53 +1543,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/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

From 8f978d95c20bf449299746031cb42a13f1bb0636 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 10 Oct 2017 09:39:58 +0200
Subject: [PATCH 14/17] Update for OCaml 4.04

---
 ocaml/Block.ml              |  4 ++--
 ocaml/Default.ml            |  2 +-
 ocaml/Input.ml              |  2 +-
 ocaml/Launcher.ml           |  2 +-
 ocaml/Md5.ml                |  2 +-
 ocaml/Message.ml            |  2 +-
 ocaml/Qmcchem_config.ml     |  2 +-
 ocaml/Qmcchem_dataserver.ml |  4 ++--
 ocaml/Qmcchem_debug.ml      |  4 ++--
 ocaml/Qmcchem_edit.ml       |  2 +-
 ocaml/Qmcchem_forwarder.ml  | 12 ++++++------
 ocaml/Qmcchem_info.ml       |  6 +++---
 ocaml/Qmcchem_md5.ml        |  2 +-
 ocaml/Qmcchem_result.ml     |  4 ++--
 ocaml/Qmcchem_run.ml        | 18 +++++++++---------
 ocaml/Qmcchem_stop.ml       |  2 +-
 ocaml/Qputils.ml            |  2 +-
 ocaml/Random_variable.ml    |  4 ++--
 ocaml/Sample.ml             |  4 ++--
 ocaml/Sample.mli            |  4 +++-
 ocaml/Scheduler.ml          |  2 +-
 ocaml/Watchdog.ml           |  6 +++---
 ocaml/build.ninja           |  4 ++--
 ocaml/ninja_ocaml.py        |  2 +-
 ocaml/qmcchem.ml            |  2 +-
 ocaml/qptypes_generator.ml  | 12 ++++++------
 26 files changed, 57 insertions(+), 55 deletions(-)

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..6386771 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 ;
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 ();

From cf003c3e08381e89351c18e52684696bebd539f5 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 10 Oct 2017 10:26:23 +0200
Subject: [PATCH 15/17] DGEMV

---
 scripts/create_properties_ezfio.py | 4 ++--
 src/det.irp.f                      | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

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/det.irp.f b/src/det.irp.f
index 17410f0..92ae46e 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -1543,7 +1543,7 @@ END_PROVIDER
   DaC = 0.d0
   CDb = 0.d0
 
-  if (det_num < ishft(det_alpha_num*det_beta_num,2)) then
+  if (det_num < ishft(det_alpha_num*det_beta_num,-2)) then
 
     det_num4 = iand(det_num,not(3))
     !DIR$ VECTOR ALIGNED

From 074bfc17058ed6b0d077a58015e067d38d53c4fd Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 12 Oct 2017 15:51:15 +0200
Subject: [PATCH 16/17] Fixed duration of a block

---
 src/SAMPLING/dmc_step.irp.f  |  2 +-
 src/SAMPLING/fkmc_step.irp.f |  2 +-
 src/SAMPLING/pdmc_step.irp.f |  2 +-
 src/SAMPLING/srmc_step.irp.f |  6 +++---
 src/det.irp.f                | 26 ++++++++------------------
 src/simulation.irp.f         |  2 +-
 6 files changed, 15 insertions(+), 25 deletions(-)

diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f
index 7044080..27397bc 100644
--- a/src/SAMPLING/dmc_step.irp.f
+++ b/src/SAMPLING/dmc_step.irp.f
@@ -248,7 +248,7 @@ END_SHELL
   if (cpu1 < cpu0) then
     cpu1 = cpu1+cpu0
   endif
-  loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time
+  loop = dble(cpu1-cpu0)/dble(count_rate) < block_time
   if (cpu1-cpu2 > count_rate) then
     integer                        :: do_run
     call get_running(do_run)
diff --git a/src/SAMPLING/fkmc_step.irp.f b/src/SAMPLING/fkmc_step.irp.f
index 3f57d20..c7e0e95 100644
--- a/src/SAMPLING/fkmc_step.irp.f
+++ b/src/SAMPLING/fkmc_step.irp.f
@@ -319,7 +319,7 @@ END_SHELL
   if (cpu1 < cpu0) then
     cpu1 = cpu1+cpu0
   endif
-  loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time
+  loop = dble(cpu1-cpu0)/dble(count_rate) < block_time
   if (cpu1-cpu2 > count_rate) then
     integer                        :: do_run
     call get_running(do_run)
diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f
index e4c3e70..484d7d8 100644
--- a/src/SAMPLING/pdmc_step.irp.f
+++ b/src/SAMPLING/pdmc_step.irp.f
@@ -280,7 +280,7 @@ END_SHELL
   if (cpu1 < cpu0) then
     cpu1 = cpu1+cpu0
   endif
-  loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time
+  loop = dble(cpu1-cpu0)/dble(count_rate) < block_time
   if (cpu1-cpu2 > count_rate) then
     integer                        :: do_run
     call get_running(do_run)
diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f
index b34692f..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,7 +320,7 @@ END_SHELL
   if (cpu1 < cpu0) then
     cpu1 = cpu1+cpu0
   endif
-  loop = dble(cpu1-cpu0)*dble(walk_num)/dble(count_rate) < block_time
+  loop = dble(cpu1-cpu0)/dble(count_rate) < block_time
   if (cpu1-cpu2 > count_rate) then
     integer                        :: do_run
     call get_running(do_run)
diff --git a/src/det.irp.f b/src/det.irp.f
index 92ae46e..9b1f6d0 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -1408,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
   
@@ -1460,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
   
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
 

From 7883be32c1d655c5c2cd369aa7ad0bc6954da2c2 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 12 Oct 2017 17:59:08 +0200
Subject: [PATCH 17/17] Removed ZMQ unbind

---
 ocaml/Qmcchem_dataserver.ml | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/ocaml/Qmcchem_dataserver.ml b/ocaml/Qmcchem_dataserver.ml
index 6386771..3067996 100644
--- a/ocaml/Qmcchem_dataserver.ml
+++ b/ocaml/Qmcchem_dataserver.ml
@@ -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