diff --git a/plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f b/plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f
index a25ff56d..e24aeb11 100644
--- a/plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f
+++ b/plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f
@@ -38,7 +38,8 @@ subroutine run_wf
   integer, external              :: zmq_get_dvector, zmq_get_N_det_generators 
   integer, external              :: zmq_get8_dvector
   integer, external              :: zmq_get_ivector
-  integer, external              :: zmq_get_psi, zmq_get_N_det_selectors
+  integer, external              :: zmq_get_psi, zmq_get_N_det_selectors, zmq_get_psi_bilinear
+  integer, external              :: zmq_get_psi_notouch
   integer, external              :: zmq_get_N_states_diag
 
   zmq_context = f77_zmq_ctx_new ()
@@ -131,8 +132,8 @@ subroutine run_wf
       ! --------
 
       call wall_time(t0)
-      if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
       if (zmq_get_N_states_diag(zmq_to_qp_run_socket,1) == -1) cycle
+      if (zmq_get_psi_bilinear(zmq_to_qp_run_socket,1) == -1) cycle
       if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states_diag) == -1) cycle
 
       call wall_time(t1)
diff --git a/src/Davidson/davidson_parallel.irp.f b/src/Davidson/davidson_parallel.irp.f
index 8d580083..5d150bb3 100644
--- a/src/Davidson/davidson_parallel.irp.f
+++ b/src/Davidson/davidson_parallel.irp.f
@@ -84,6 +84,11 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, N_st, sze,
   integer, external :: zmq_get_dvector
   integer, external :: zmq_get_dmatrix
 
+  PROVIDE psi_det_beta_unique psi_bilinear_matrix_order_transp_reverse psi_det_alpha_unique 
+  PROVIDE psi_bilinear_matrix_transp_values psi_bilinear_matrix_values psi_bilinear_matrix_columns_loc
+  PROVIDE ref_bitmask_energy nproc
+  PROVIDE mpi_initialized
+
   allocate(u_t(N_st,N_det))
   allocate (energy(N_st))
 
@@ -96,11 +101,11 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, N_st, sze,
     ni = 8388608
     nj = int(size(u_t,kind=8)/8388608_8,4) + 1
   endif
-  if (zmq_get_dmatrix(zmq_to_qp_run_socket, worker_id, 'u_t', u_t, ni, nj, size(u_t,kind=8)) == -1) then
-    print *,  irp_here, ': Unable to get u_t'
-    deallocate(u_t,energy)
-    return
-  endif
+
+  do while (zmq_get_dmatrix(zmq_to_qp_run_socket, worker_id, 'u_t', u_t, ni, nj, size(u_t,kind=8)) == -1) 
+    call sleep(1)
+    print *,  irp_here, ': waiting for u_t...'
+  enddo
 
   if (zmq_get_dvector(zmq_to_qp_run_socket, worker_id, 'energy', energy, size(energy)) == -1) then
     print *,  irp_here, ': Unable to get energy'
@@ -293,17 +298,71 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
   integer                        :: ithread
   double precision, allocatable  :: u_t(:,:)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
-  
+  integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
   PROVIDE psi_det_beta_unique psi_bilinear_matrix_order_transp_reverse psi_det_alpha_unique 
   PROVIDE psi_bilinear_matrix_transp_values psi_bilinear_matrix_values psi_bilinear_matrix_columns_loc
   PROVIDE ref_bitmask_energy nproc
   PROVIDE mpi_initialized
 
+  call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'davidson')
+  
+  if (zmq_put_N_states_diag(zmq_to_qp_run_socket, 1) == -1) then
+    stop 'Unable to put N_states_diag on ZMQ server'
+  endif
+  if (zmq_put_psi_bilinear(zmq_to_qp_run_socket,1) == -1) then
+    stop 'Unable to put psi on ZMQ server'
+  endif
+  energy = 0.d0
+  if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',energy,size(energy)) == -1) then
+    stop 'Unable to put energy on ZMQ server'
+  endif
+
+
+  ! Create tasks
+  ! ============
+
+  integer :: istep, imin, imax, ishift, ipos
+  integer, external :: add_task_to_taskserver
+  integer, parameter :: tasksize=10000
+  character*(100000) :: task
+  istep=1
+  ishift=0
+  imin=1
+
+
+  ipos=1
+  do imin=1,N_det,10000
+    imax = min(N_det,imin-1+tasksize)
+    do ishift=0,istep-1
+      write(task(ipos:ipos+50),'(4(I11,1X),1X,1A)') imin, imax, ishift, istep, '|'
+      ipos = ipos+50
+      if (ipos > 100000-50) then
+        if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
+          stop 'Unable to add task'
+        endif
+        ipos=1
+      endif
+    enddo
+  enddo
+
+  if (ipos > 1) then
+    if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
+      stop 'Unable to add task'
+    endif
+    ipos=1
+  endif
+    
+  integer, external :: zmq_set_running
+  if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
+    print *,  irp_here, ': Failed in zmq_set_running'
+  endif
+
 
   allocate(u_t(N_st,N_det))
   do k=1,N_st
     call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
   enddo
+
   call dtranspose(                                                   &
       u_0,                                                           &
       size(u_0, 1),                                                  &
@@ -311,33 +370,17 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
       size(u_t, 1),                                                  &
       N_det, N_st)
 
-
-  integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
   
   ASSERT (N_st == N_states_diag)
   ASSERT (sze >= N_det) 
 
-  call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'davidson')
-  
-  character*(512) :: task
   integer :: rc, ni, nj
   integer*8 :: rc8
   double precision :: energy(N_st)
 
-  integer, external :: zmq_put_dvector, zmq_put_psi, zmq_put_N_states_diag
+  integer, external :: zmq_put_dvector, zmq_put_psi, zmq_put_N_states_diag, zmq_put_psi_bilinear
   integer, external :: zmq_put_dmatrix
 
-  energy = 0.d0
-
-  if (zmq_put_N_states_diag(zmq_to_qp_run_socket, 1) == -1) then
-    stop 'Unable to put N_states_diag on ZMQ server'
-  endif
-  if (zmq_put_psi(zmq_to_qp_run_socket,1) == -1) then
-    stop 'Unable to put psi on ZMQ server'
-  endif
-  if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',energy,size(energy)) == -1) then
-    stop 'Unable to put energy on ZMQ server'
-  endif
   if (size(u_t) < 8388608) then
     ni = size(u_t)
     nj = 1
@@ -354,50 +397,9 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
   deallocate(u_t)
 
 
-  ! Create tasks
-  ! ============
-
-  integer :: istep, imin, imax, ishift
-  double precision :: w, max_workload, N_det_inv
-  integer, external :: add_task_to_taskserver
-  w = 0.d0
-  istep=1
-  ishift=0
-  imin=1
-  N_det_inv = 1.d0/dble(N_det)
-  max_workload = 10000.d0
-  do imax=1,N_det
-    w = w + 1.d0
-    if (w > max_workload) then
-      do ishift=0,istep-1
-        write(task,'(4(I9,1X),1A)') imin, imax, ishift, istep, '|'
-        if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then
-          stop 'Unable to add task'
-        endif
-      enddo
-      imin = imax+1
-      w = 0.d0
-    endif
-  enddo
-  if (w > 0.d0) then
-    imax = N_det
-    do ishift=0,istep-1
-      write(task,'(4(I9,1X),1A)') imin, imax, ishift, istep, '|'
-      if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then
-        stop 'Unable to add task'
-      endif
-    enddo
-  endif
-    
-
   v_0 = 0.d0
   s_0 = 0.d0
 
-  integer, external :: zmq_set_running
-  if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
-    print *,  irp_here, ': Failed in zmq_set_running'
-  endif
-
   call omp_set_nested(.True.)
   !$OMP PARALLEL NUM_THREADS(2) PRIVATE(ithread)
   ithread = omp_get_thread_num()
diff --git a/src/DavidsonUndressed/print_energy.irp.f b/src/DavidsonUndressed/print_energy.irp.f
index d694cb6c..94165560 100644
--- a/src/DavidsonUndressed/print_energy.irp.f
+++ b/src/DavidsonUndressed/print_energy.irp.f
@@ -2,7 +2,7 @@ program print_energy
  implicit none
  read_wf = .true.
  touch read_wf
- provide mo_bielec_integrals_in_map
+ provide mo_bielec_integrals_in_map psi_coef psi_det psi_bilinear_matrix_transp_values
  double precision :: time1, time0
  call wall_time(time0)
  call routine
diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f
index 9a39212c..29d0eb30 100644
--- a/src/Determinants/spindeterminants.irp.f
+++ b/src/Determinants/spindeterminants.irp.f
@@ -519,7 +519,7 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_
     do k=1,N_det
       psi_bilinear_matrix_transp_values (k,l) = psi_bilinear_matrix_values (k,l)
     enddo
-    !$OMP ENDDO
+    !$OMP ENDDO NOWAIT
   enddo
   !$OMP DO 
   do k=1,N_det
diff --git a/src/Determinants/zmq.irp.f b/src/Determinants/zmq.irp.f
index 97af6210..33350974 100644
--- a/src/Determinants/zmq.irp.f
+++ b/src/Determinants/zmq.irp.f
@@ -40,6 +40,242 @@ end
 
 
 
+integer function zmq_get_psi_notouch(zmq_to_qp_run_socket, worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get the wave function from the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+
+  integer, external              :: zmq_get_N_states
+  integer, external              :: zmq_get_N_det
+  integer, external              :: zmq_get_psi_det_size
+  integer*8, external            :: zmq_get_psi_det
+  integer*8, external            :: zmq_get_psi_coef
+
+  zmq_get_psi_notouch = 0
+
+  if (zmq_get_N_states(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_get_psi_notouch = -1
+    return
+  endif
+  if (zmq_get_N_det(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_get_psi_notouch = -1
+    return
+  endif
+  if (zmq_get_psi_det_size(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_get_psi_notouch = -1
+    return
+  endif
+  
+  if (size(psi_det,kind=8) /= N_int*2_8*psi_det_size*bit_kind) then
+    deallocate(psi_det)
+    allocate(psi_det(N_int,2,psi_det_size))
+  endif
+
+  if (size(psi_coef,kind=8) /= psi_det_size*N_states) then
+    deallocate(psi_coef)
+    allocate(psi_coef(psi_det_size,N_states))
+  endif
+
+  if (zmq_get_psi_det(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_notouch = -1
+    return
+  endif
+  if (zmq_get_psi_coef(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_notouch = -1
+    return
+  endif
+
+end
+
+
+integer function zmq_get_psi(zmq_to_qp_run_socket, worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get the wave function from the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer, external :: zmq_get_psi_notouch
+  zmq_get_psi = zmq_get_psi_notouch(zmq_to_qp_run_socket, worker_id)
+  SOFT_TOUCH psi_det psi_coef psi_det_size N_det N_states
+
+end
+
+
+
+
+
+integer function zmq_put_psi_bilinear(zmq_to_qp_run_socket,worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Put the wave function on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  character*(256)                :: msg
+
+
+  zmq_put_psi_bilinear = 0
+
+  integer, external            :: zmq_put_psi
+  if (zmq_put_psi(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_put_psi_bilinear_matrix_columns
+  if (zmq_put_psi_bilinear_matrix_rows(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_put_psi_bilinear_matrix_rows
+  if (zmq_put_psi_bilinear_matrix_columns(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_put_psi_bilinear_matrix_order
+  if (zmq_put_psi_bilinear_matrix_order(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_put_psi_bilinear_matrix_values
+  if (zmq_put_psi_bilinear_matrix_values(zmq_to_qp_run_socket, worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_put_N_det_alpha_unique
+  if (zmq_put_N_det_alpha_unique(zmq_to_qp_run_socket,worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_put_N_det_beta_unique
+  if (zmq_put_N_det_beta_unique(zmq_to_qp_run_socket,worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_put_psi_det_alpha_unique
+  if (zmq_put_psi_det_alpha_unique(zmq_to_qp_run_socket,worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_put_psi_det_beta_unique
+  if (zmq_put_psi_det_beta_unique(zmq_to_qp_run_socket,worker_id) == -1) then
+    zmq_put_psi_bilinear = -1
+    return
+  endif
+
+end
+
+
+integer function zmq_get_psi_bilinear(zmq_to_qp_run_socket, worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get the wave function from the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+
+  integer, external            :: zmq_get_psi_notouch
+  if (zmq_get_psi_notouch(zmq_to_qp_run_socket,1) == -1) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+
+  zmq_get_psi_bilinear= 0
+
+  if (size(psi_bilinear_matrix_values,kind=8) /= N_det*N_states) then
+    deallocate(psi_bilinear_matrix_values)
+    allocate(psi_bilinear_matrix_values(N_det,N_states))
+  endif
+
+  if (size(psi_bilinear_matrix_rows,kind=8) /= N_det) then
+    deallocate(psi_bilinear_matrix_rows)
+    allocate(psi_bilinear_matrix_rows(N_det))
+  endif
+
+  if (size(psi_bilinear_matrix_columns,kind=8) /= N_det) then
+    deallocate(psi_bilinear_matrix_columns)
+    allocate(psi_bilinear_matrix_columns(N_det))
+  endif
+
+  if (size(psi_bilinear_matrix_order,kind=8) /= N_det) then
+    deallocate(psi_bilinear_matrix_order)
+    allocate(psi_bilinear_matrix_order(N_det))
+  endif
+
+  integer*8, external            :: zmq_get_psi_bilinear_matrix_values
+  if (zmq_get_psi_bilinear_matrix_values(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_get_psi_bilinear_matrix_rows
+  if (zmq_get_psi_bilinear_matrix_rows(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_get_psi_bilinear_matrix_columns
+  if (zmq_get_psi_bilinear_matrix_columns(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_get_psi_bilinear_matrix_order
+  if (zmq_get_psi_bilinear_matrix_order(zmq_to_qp_run_socket, worker_id) == -1_8) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+
+
+  integer*8, external            :: zmq_get_N_det_alpha_unique
+  if (zmq_get_N_det_alpha_unique(zmq_to_qp_run_socket,worker_id) == -1) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_get_N_det_beta_unique
+  if (zmq_get_N_det_beta_unique(zmq_to_qp_run_socket,worker_id) == -1) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_get_psi_det_alpha_unique
+  if (zmq_get_psi_det_alpha_unique(zmq_to_qp_run_socket,worker_id) == -1) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+
+  integer*8, external            :: zmq_get_psi_det_beta_unique
+  if (zmq_get_psi_det_beta_unique(zmq_to_qp_run_socket,worker_id) == -1) then
+    zmq_get_psi_bilinear = -1
+    return
+  endif
+
+  SOFT_TOUCH psi_bilinear_matrix_values psi_bilinear_matrix_rows psi_bilinear_matrix_columns psi_bilinear_matrix_order psi_det psi_coef psi_det_size N_det N_states psi_det_beta_unique psi_det_alpha_unique N_det_beta_unique N_det_alpha_unique
+
+end
+
+
+
+
+
+
+
 BEGIN_TEMPLATE 
 
 integer function zmq_put_$X(zmq_to_qp_run_socket,worker_id)
@@ -136,15 +372,20 @@ SUBST [ X ]
 
 N_states ;;
 N_det ;;
+N_det_alpha_unique ;;
+N_det_beta_unique ;;
 psi_det_size ;;
 
 END_TEMPLATE
 
-integer*8 function zmq_put_psi_det(zmq_to_qp_run_socket,worker_id)
+
+BEGIN_TEMPLATE
+
+integer*8 function zmq_put_$X(zmq_to_qp_run_socket,worker_id)
   use f77_zmq
   implicit none
   BEGIN_DOC
-! Put psi_det on the qp_run scheduler
+! Put $X on the qp_run scheduler
   END_DOC
   integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
   integer, intent(in)            :: worker_id
@@ -154,49 +395,22 @@ integer*8 function zmq_put_psi_det(zmq_to_qp_run_socket,worker_id)
   integer*8 :: zmq_put_i8matrix
   integer   :: ni, nj
 
-  if (size(psi_det,kind=8) <= 8388608_8) then
-    ni = size(psi_det,kind=4)
+  if (size($X,kind=8) <= 8388608_8) then
+    ni = size($X,kind=4)
     nj = 1
   else
     ni = 8388608_8
-    nj = int(size(psi_det,kind=8)/8388608_8,4) + 1
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
   endif
-  rc8 = zmq_put_i8matrix(zmq_to_qp_run_socket, 1, 'psi_det', psi_det, ni, nj, size(psi_det,kind=8))
-  zmq_put_psi_det = rc8
+  rc8 = zmq_put_i8matrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8))
+  zmq_put_$X = rc8
 end
 
-integer*8 function zmq_put_psi_coef(zmq_to_qp_run_socket,worker_id)
+integer*8 function zmq_get_$X(zmq_to_qp_run_socket,worker_id)
   use f77_zmq
   implicit none
   BEGIN_DOC
-! Put psi_coef on the qp_run scheduler
-  END_DOC
-  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
-  integer, intent(in)            :: worker_id
-  integer*8                      :: rc8
-  character*(256)                :: msg
-
-  zmq_put_psi_coef = 0
-
-  integer*8 :: zmq_put_dmatrix
-  integer   :: ni, nj
-
-  if (size(psi_coef,kind=8) <= 8388608_8) then
-    ni = size(psi_coef,kind=4)
-    nj = 1
-  else
-    ni = 8388608
-    nj = int(size(psi_coef,kind=8)/8388608_8,4) + 1
-  endif
-  rc8 = zmq_put_dmatrix(zmq_to_qp_run_socket, 1, 'psi_coef', psi_coef, ni, nj, size(psi_coef,kind=8) )
-  zmq_put_psi_coef = rc8
-end
-
-integer*8 function zmq_get_psi_det(zmq_to_qp_run_socket,worker_id)
-  use f77_zmq
-  implicit none
-  BEGIN_DOC
-! Get psi_det on the qp_run scheduler
+! Get $X on the qp_run scheduler
   END_DOC
   integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
   integer, intent(in)            :: worker_id
@@ -206,99 +420,150 @@ integer*8 function zmq_get_psi_det(zmq_to_qp_run_socket,worker_id)
   integer*8 :: zmq_get_i8matrix
   integer   :: ni, nj
 
-  if (size(psi_det,kind=8) <= 8388608_8) then
-    ni = size(psi_det,kind=4)
+  if (size($X,kind=8) <= 8388608_8) then
+    ni = size($X,kind=4)
     nj = 1
   else
     ni = 8388608
-    nj = int(size(psi_det,kind=8)/8388608_8,4) + 1
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
   endif
-  rc8 = zmq_get_i8matrix(zmq_to_qp_run_socket, 1, 'psi_det', psi_det, ni, nj, size(psi_det,kind=8))
-  zmq_get_psi_det = rc8
+  rc8 = zmq_get_i8matrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8))
+  zmq_get_$X = rc8
 end
 
-integer*8 function zmq_get_psi_coef(zmq_to_qp_run_socket,worker_id)
+SUBST [ X ]
+
+psi_det ;;
+psi_det_alpha_unique ;;
+psi_det_beta_unique ;;
+
+END_TEMPLATE
+
+BEGIN_TEMPLATE
+
+integer*8 function zmq_put_$X(zmq_to_qp_run_socket,worker_id)
   use f77_zmq
   implicit none
   BEGIN_DOC
-! get psi_coef on the qp_run scheduler
+! Put $X on the qp_run scheduler
   END_DOC
   integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
   integer, intent(in)            :: worker_id
   integer*8                      :: rc8
   character*(256)                :: msg
 
-  zmq_get_psi_coef = 0_8
+  integer*8 :: zmq_put_imatrix
+  integer   :: ni, nj
+
+  if (size($X,kind=8) <= 8388608_8) then
+    ni = size($X,kind=4)
+    nj = 1
+  else
+    ni = 8388608_8
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
+  endif
+  rc8 = zmq_put_imatrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8))
+  zmq_put_$X = rc8
+end
+
+integer*8 function zmq_get_$X(zmq_to_qp_run_socket,worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get $X on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer*8                      :: rc8
+  character*(256)                :: msg
+
+  integer*8 :: zmq_get_imatrix
+  integer   :: ni, nj
+
+  if (size($X,kind=8) <= 8388608_8) then
+    ni = size($X,kind=4)
+    nj = 1
+  else
+    ni = 8388608
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
+  endif
+  rc8 = zmq_get_imatrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8))
+  zmq_get_$X = rc8
+end
+
+SUBST [ X ]
+
+psi_bilinear_matrix_rows ;;
+psi_bilinear_matrix_columns ;;
+psi_bilinear_matrix_order ;;
+
+END_TEMPLATE
+
+
+BEGIN_TEMPLATE
+
+integer*8 function zmq_put_$X(zmq_to_qp_run_socket,worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Put $X on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer*8                      :: rc8
+  character*(256)                :: msg
+
+  zmq_put_$X = 0
+
+  integer*8 :: zmq_put_dmatrix
+  integer   :: ni, nj
+
+  if (size($X,kind=8) <= 8388608_8) then
+    ni = size($X,kind=4)
+    nj = 1
+  else
+    ni = 8388608
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
+  endif
+  rc8 = zmq_put_dmatrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8) )
+  zmq_put_$X = rc8
+end
+
+integer*8 function zmq_get_$X(zmq_to_qp_run_socket,worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! get $X on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer*8                      :: rc8
+  character*(256)                :: msg
+
+  zmq_get_$X = 0_8
 
   integer*8 :: zmq_get_dmatrix
   integer   :: ni, nj
 
-  if (size(psi_coef,kind=8) <= 8388608_8) then
-    ni = size(psi_coef,kind=4)
+  if (size($X,kind=8) <= 8388608_8) then
+    ni = size($X,kind=4)
     nj = 1
   else
     ni = 8388608
-    nj = int(size(psi_coef,kind=8)/8388608_8,4) + 1
+    nj = int(size($X,kind=8)/8388608_8,4) + 1
   endif
-  rc8 = zmq_get_dmatrix(zmq_to_qp_run_socket, 1, 'psi_coef', psi_coef, ni, nj, size(psi_coef,kind=8) )
-  zmq_get_psi_coef = rc8
+  rc8 = zmq_get_dmatrix(zmq_to_qp_run_socket, 1, '$X', $X, ni, nj, size($X,kind=8) )
+  zmq_get_$X = rc8
 end
 
+SUBST [ X ]
+
+psi_coef ;;
+psi_bilinear_matrix_values ;;
+
+END_TEMPLATE
+
+
 !---------------------------------------------------------------------------
 
 
-integer function zmq_get_psi(zmq_to_qp_run_socket, worker_id)
-  use f77_zmq
-  implicit none
-  BEGIN_DOC
-! Get the wave function from the qp_run scheduler
-  END_DOC
-  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
-  integer, intent(in)            :: worker_id
-
-  integer, external              :: zmq_get_N_states
-  integer, external              :: zmq_get_N_det
-  integer, external              :: zmq_get_psi_det_size
-  integer*8, external            :: zmq_get_psi_det
-  integer*8, external            :: zmq_get_psi_coef
-
-  zmq_get_psi = 0
-
-  if (zmq_get_N_states(zmq_to_qp_run_socket, worker_id) == -1) then
-    zmq_get_psi = -1
-    return
-  endif
-  if (zmq_get_N_det(zmq_to_qp_run_socket, worker_id) == -1) then
-    zmq_get_psi = -1
-    return
-  endif
-  if (zmq_get_psi_det_size(zmq_to_qp_run_socket, worker_id) == -1) then
-    zmq_get_psi = -1
-    return
-  endif
-  
-  if (size(psi_det,kind=8) /= N_int*2_8*psi_det_size*bit_kind) then
-    deallocate(psi_det)
-    allocate(psi_det(N_int,2,psi_det_size))
-  endif
-
-  if (size(psi_coef,kind=8) /= psi_det_size*N_states) then
-    deallocate(psi_coef)
-    allocate(psi_coef(psi_det_size,N_states))
-  endif
-
-  if (zmq_get_psi_det(zmq_to_qp_run_socket, worker_id) == -1_8) then
-    zmq_get_psi = -1
-    return
-  endif
-  if (zmq_get_psi_coef(zmq_to_qp_run_socket, worker_id) == -1_8) then
-    zmq_get_psi = -1
-    return
-  endif
-  SOFT_TOUCH psi_det psi_coef psi_det_size N_det N_states
-
-end
-
-
-
-
diff --git a/src/Utils/sort.irp.f b/src/Utils/sort.irp.f
index a9594d6c..085d3d35 100644
--- a/src/Utils/sort.irp.f
+++ b/src/Utils/sort.irp.f
@@ -456,16 +456,16 @@ BEGIN_TEMPLATE
         iorder(i) = iorder1(1_$int_type+i1-i)
       enddo
     endif
-    deallocate(x1,iorder1,stat=err)
-    if (err /= 0) then
-      print *,  irp_here, ': Unable to deallocate arrays x1, iorder1'
-      stop
-    endif
     
     if (i2>1_$int_type) then
       call $Xradix_sort$big(x(i1+1_$int_type),iorder(i1+1_$int_type),i2,-2)
     endif
     
+    deallocate(x1,iorder1,stat=err)
+    if (err /= 0) then
+      print *,  irp_here, ': Unable to deallocate arrays x1, iorder1'
+      stop
+    endif
     return
 
   else if (iradix == -2) then ! Positive
@@ -526,13 +526,23 @@ BEGIN_TEMPLATE
     endif
     
     
+    !$OMP PARALLEL DEFAULT(SHARED) if (isize > 1000000)
+    !$OMP SINGLE
     if (i3>1_$int_type) then
+      !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(i3 > 1000000)
       call $Xradix_sort$big(x,iorder,i3,iradix_new-1)
+      !$OMP END TASK
     endif
     
     if (isize-i3>1_$int_type) then
+      !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(isize-i3 > 1000000)
       call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1)
+      !$OMP END TASK
     endif
+
+    !$OMP TASKWAIT
+    !$OMP END SINGLE
+    !$OMP END PARALLEL
     
     return
   endif
@@ -588,11 +598,16 @@ BEGIN_TEMPLATE
   
   
   if (i1>1_$int_type) then
+    !$OMP TASK FIRSTPRIVATE(i0,iradix,i1) SHARED(x,iorder) if(i1 >1000000)
     call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1)
+    !$OMP END TASK
   endif
   if (i0>1) then
+    !$OMP TASK FIRSTPRIVATE(i0,iradix) SHARED(x,iorder) if(i0 >1000000)
     call $Xradix_sort$big(x,iorder,i0,iradix-1)
+    !$OMP END TASK
   endif
+  !$OMP TASKWAIT
   
  end
 
diff --git a/src/ZMQ/put_get.irp.f b/src/ZMQ/put_get.irp.f
index 4fb4ecea..95a4ef4a 100644
--- a/src/ZMQ/put_get.irp.f
+++ b/src/ZMQ/put_get.irp.f
@@ -608,13 +608,13 @@ integer function zmq_get_int(zmq_to_qp_run_socket, worker_id, name, x)
     include 'mpif.h'
     call MPI_BCAST (zmq_get_int, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
     if (ierr /= MPI_SUCCESS) then
-      print *,  irp_here//': Unable to broadcast zmq_get_i8matrix'
+      print *,  irp_here//': Unable to broadcast zmq_get_int'
       stop -1
     endif
     call MPI_BARRIER(MPI_COMM_WORLD,ierr)
     call MPI_BCAST (x, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
     if (ierr /= MPI_SUCCESS) then
-      print *,  irp_here//': Unable to broadcast zmq_get_i8matrix'
+      print *,  irp_here//': Unable to broadcast zmq_get_int'
       stop -1
     endif
   IRP_ENDIF
@@ -723,7 +723,7 @@ integer function zmq_get_i8matrix(zmq_to_qp_run_socket, worker_id, name, x, size
   integer, intent(in)            :: size_x1, size_x2
   integer*8, intent(in)          :: sze
   character*(*), intent(in)      :: name
-  double precision, intent(out)  :: x(size_x1,size_x2)
+  integer*8, intent(out)         :: x(size_x1,size_x2)
   integer*8                      :: rc, ni
   integer*8                      :: j
   character*(256)                :: msg
@@ -777,7 +777,7 @@ integer function zmq_get_i8matrix(zmq_to_qp_run_socket, worker_id, name, x, size
       stop -1
     endif
     call MPI_BARRIER(MPI_COMM_WORLD,ierr)
-    call broadcast_chunks_double(x, sze)
+    call broadcast_chunks_integer8(x, sze)
   IRP_ENDIF
 
 end
@@ -786,3 +786,124 @@ end
 
 
 
+integer function zmq_put_imatrix(zmq_to_qp_run_socket, worker_id, name, x, size_x1, size_x2, sze)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Put a float vector on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  character*(*)                  :: name
+  integer, intent(in)            :: size_x1, size_x2
+  integer*8, intent(in)          :: sze
+  integer, intent(in)            :: x(size_x1, size_x2)
+  integer*8                      :: rc, ni
+  integer*8                      :: j
+  character*(256)                :: msg
+
+  zmq_put_imatrix = 0
+
+  ni = size_x1
+  do j=1,size_x2
+    if (j == size_x2) then
+      ni = sze - (j-1_8)*int(size_x1,8)
+    endif
+    write(msg,'(A,1X,I8,1X,A,I8.8)') 'put_data '//trim(zmq_state), worker_id, trim(name), j
+    rc = f77_zmq_send(zmq_to_qp_run_socket,trim(msg),len(trim(msg)),ZMQ_SNDMORE)
+    if (rc /= len(trim(msg))) then
+      zmq_put_imatrix = -1
+      print *,  irp_here, 'Failed in put_data', rc, j
+      return 
+    endif
+
+    rc = f77_zmq_send8(zmq_to_qp_run_socket,x(1,j),ni*4_8,0)
+    if (rc /= ni*4_8) then
+      print *,  irp_here, 'Failed in send ', rc, j
+      zmq_put_imatrix = -1
+      return
+    endif
+
+    rc = f77_zmq_recv(zmq_to_qp_run_socket,msg,len(msg),0)
+    if (msg(1:rc) /= 'put_data_reply ok') then
+      print *,  irp_here, 'Failed in recv ',  rc, j
+      zmq_put_imatrix = -1
+      return
+    endif
+  enddo
+
+end
+
+
+integer function zmq_get_imatrix(zmq_to_qp_run_socket, worker_id, name, x, size_x1, size_x2, sze)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get a float vector from the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer, intent(in)            :: size_x1, size_x2
+  integer*8, intent(in)          :: sze
+  character*(*), intent(in)      :: name
+  integer, intent(out)           :: x(size_x1,size_x2)
+  integer*8                      :: rc, ni
+  integer*8                      :: j
+  character*(256)                :: msg
+
+  PROVIDE zmq_state
+  ! Success
+  zmq_get_imatrix = 0
+
+  if (mpi_master) then
+    ni = size_x1
+    do j=1, size_x2
+      if (j == size_x2) then
+        ni = sze - (j-1)*size_x1
+      endif
+      write(msg,'(A,1X,I8,1X,A,I8.8)') 'get_data '//trim(zmq_state), worker_id, trim(name),j
+      rc = f77_zmq_send(zmq_to_qp_run_socket,trim(msg),len(trim(msg)),0)
+      if (rc /= len(trim(msg))) then
+        zmq_get_imatrix = -1
+        print *,  irp_here, 'rc /= len(trim(msg))', rc, len(trim(msg))
+        go to 10
+      endif
+
+      rc = f77_zmq_recv(zmq_to_qp_run_socket,msg,len(msg),0)
+      if (msg(1:14) /= 'get_data_reply') then
+        print *,  irp_here, 'msg(1:14) /= get_data_reply', msg(1:14)
+        zmq_get_imatrix = -1
+        go to 10
+      endif
+
+      rc = f77_zmq_recv8(zmq_to_qp_run_socket,x(1,j),ni*4_8,0)
+      if (rc /= ni*4_8) then
+        print *,  irp_here, 'rc /= ni*8', rc, ni*4_8
+        zmq_get_imatrix = -1
+        go to 10
+      endif
+    enddo
+  endif
+
+  10 continue
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+    integer :: ierr
+    include 'mpif.h'
+    call MPI_BCAST (zmq_get_imatrix, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      print *,  irp_here//': Unable to broadcast zmq_get_imatrix'
+      stop -1
+    endif
+    call MPI_BARRIER(MPI_COMM_WORLD,ierr)
+    call broadcast_chunks_integer(x, sze)
+  IRP_ENDIF
+
+end
+
+
+