Merge branch 'dev' of github.com:QuantumPackage/qp2 into dev

src/ao_two_e_ints/integrals_in_map_slave.irp.f

@@ -116,7 +116,7 @@ subroutine ao_two_e_integrals_in_map_slave(thread,iproc)
       exit
     endif
     if (task_id == 0) exit
-    read(task,*) j, l
+    call sscanf_dd(task, j, l)
     integer, external :: task_done_to_taskserver
     call compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value)
     if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) == -1) then

src/cipsi/pt2_stoch_routines.irp.f

@@ -187,7 +187,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
       if (zmq_put_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) then
         stop 'Unable to put pt2_stoch_istate on ZMQ server'
       endif
-      if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) then
+      if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',(/threshold_generators/),1) == -1) then
         stop 'Unable to put threshold_generators on ZMQ server'
       endif
 

src/cipsi/run_selection_slave.irp.f

@@ -52,7 +52,7 @@ subroutine run_selection_slave(thread,iproc,energy)
       ctask = ctask - 1
     else
       integer :: i_generator, N, subset, bsize
-      read(task,*) subset, i_generator, N
+      call sscanf_ddd(task, subset, i_generator, N)
       if(buf%N == 0) then
         ! Only first time
         call create_selection_buffer(N, N*2, buf)

src/cipsi/selection.irp.f

@@ -8,12 +8,15 @@ subroutine get_mask_phase(det1, pm, Nint)
   integer(bit_kind), intent(out) :: pm(Nint,2)
   integer(bit_kind) :: tmp1, tmp2
   integer :: i
-  pm(1:Nint,1:2) = det1(1:Nint,1:2)
   tmp1 = 0_8
   tmp2 = 0_8
-  do i=1,Nint
+  select case (Nint)
-    pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 1))
+
-    pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 1))
+BEGIN_TEMPLATE
+    case ($Nint)
+      do i=1,$Nint
+        pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1))
+        pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1))
         pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2))
         pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2))
         pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4))
@@ -29,6 +32,32 @@ subroutine get_mask_phase(det1, pm, Nint)
         if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1)
         if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2)
       end do
+SUBST [ Nint ]
+1;;
+2;;
+3;;
+4;;
+END_TEMPLATE
+    case default
+      do i=1,Nint
+        pm(i,1) = ieor(det1(i,1), shiftl(det1(i,1), 1))
+        pm(i,2) = ieor(det1(i,2), shiftl(det1(i,2), 1))
+        pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 2))
+        pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 2))
+        pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 4))
+        pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 4))
+        pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 8))
+        pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 8))
+        pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 16))
+        pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 16))
+        pm(i,1) = ieor(pm(i,1), shiftl(pm(i,1), 32))
+        pm(i,2) = ieor(pm(i,2), shiftl(pm(i,2), 32))
+        pm(i,1) = ieor(pm(i,1), tmp1)
+        pm(i,2) = ieor(pm(i,2), tmp2)
+        if(iand(popcnt(det1(i,1)), 1) == 1) tmp1 = not(tmp1)
+        if(iand(popcnt(det1(i,2)), 1) == 1) tmp2 = not(tmp2)
+      end do
+  end select
 
 end subroutine
 
@@ -445,11 +474,17 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
       endif
 
       do i=1,fullinteresting(0)
-        fullminilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,fullinteresting(i))
+        do k=1,N_int
+          fullminilist(k,1,i) = psi_det_sorted(k,1,fullinteresting(i))
+          fullminilist(k,2,i) = psi_det_sorted(k,2,fullinteresting(i))
+        enddo
       enddo
 
       do i=1,interesting(0)
-        minilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,interesting(i))
+        do k=1,N_int
+          minilist(k,1,i) = psi_det_sorted(k,1,interesting(i))
+          minilist(k,2,i) = psi_det_sorted(k,2,interesting(i))
+        enddo
       enddo
 
       do s2=s1,2
@@ -668,10 +703,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 
       w = 0d0
 
-!      integer(bit_kind) :: occ(N_int,2), n
-!      call configuration_of_det(det,occ,N_int)
-!      call configuration_to_dets_size(occ,n,elec_alpha_num,N_int)
-
       e_pert = 0.d0
       coef = 0.d0
       logical :: do_diag
@@ -699,7 +730,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 
       double precision :: eigvalues(N_states+1)
       double precision :: work(1+6*(N_states+1)+2*(N_states+1)**2)
-      integer :: iwork(3+5*(N_states+1)), info, k ,n
+      integer :: iwork(3+5*(N_states+1)), info, k 
 
       if (do_diag) then
         double precision :: pt2_matrix(N_states+1,N_states+1)
@@ -765,36 +796,43 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 
           case(5)
             ! Variance selection
-!            w = w - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)
+            if (h0_type == 'CFG') then
+              w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)) & 
+                / c0_weight(istate)
+            else
               w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate))
-!            do jstate=1,N_states
+            endif
-!              if (istate == jstate) cycle
-!              w = w + dabs(alpha_h_psi*mat(jstate,p1,p2)) * s_weight(istate,jstate)
-!            enddo
 
           case(6)
-!            w = w - coef(istate) * coef(istate) * s_weight(istate,istate)
+            if (h0_type == 'CFG') then
+              w = min(w,- coef(istate) * coef(istate) * s_weight(istate,istate)) &
+                / c0_weight(istate)
+            else
               w = min(w,- coef(istate) * coef(istate) * s_weight(istate,istate))
-!            do jstate=1,N_states
+            endif
-!              if (istate == jstate) cycle
-!              w = w + dabs(coef(istate)*coef(jstate)) * s_weight(istate,jstate)
-!            enddo
 
           case default
             ! Energy selection
-!            w = w + e_pert(istate) * s_weight(istate,istate)
+            if (h0_type == 'CFG') then
+              w = min(w, e_pert(istate) * s_weight(istate,istate)) / c0_weight(istate)
+            else
               w = min(w, e_pert(istate) * s_weight(istate,istate))
-!            do jstate=1,N_states
+            endif
-!              if (istate == jstate) cycle
-!              w = w + dabs(X(istate)*X(jstate)) * s_weight(istate,jstate)
-!            enddo
 
         end select
       end do
 
 
-!      w = dble(n) * w
+      integer(bit_kind) :: occ(N_int,2), n
-
+      if (h0_type == 'CFG') then
+        do k=1,N_int
+          occ(k,1) = ieor(det(k,1),det(k,2))
+          occ(k,2) = iand(det(k,1),det(k,2))
+        enddo
+        call configuration_to_dets_size(occ,n,elec_alpha_num,N_int)
+        n = max(n,1)
+        w *= dble(n)
+      endif
 
       if(w <= buf%mini) then
         call add_to_selection_buffer(buf, det, w)

src/davidson/davidson_parallel.irp.f

@@ -438,7 +438,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
   ipos=1
   do imin=1,N_det,tasksize
     imax = min(N_det,imin-1+tasksize)
-    if (imin==1) then
+    if (imin<=N_det/2) then
       istep = 2
     else
       istep = 1
@@ -505,7 +505,9 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
     print *,  irp_here, ': Failed in zmq_set_running'
   endif
 
-  call omp_set_max_active_levels(4)
+
+  call omp_set_max_active_levels(5)
+
   !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread)
   ithread = omp_get_thread_num()
   if (ithread == 0 ) then

src/determinants/spindeterminants.irp.f

@@ -1075,19 +1075,17 @@ subroutine get_all_spin_singles_and_doubles_1(buffer, idx, spindet, size_buffer,
   integer                        :: i
   include 'utils/constants.include.F'
   integer                        :: degree
-
+  integer                        :: add_double(0:64) = (/ 0, 0, 0, 0, 1, (0, i=1,60) /)
+  integer                        :: add_single(0:64) = (/ 0, 0, 1, 0, 0, (0, i=1,60) /)
 
   n_singles = 1
   n_doubles = 1
   do i=1,size_buffer
     degree =  popcnt(  xor( spindet, buffer(i) ) )
-    if ( degree == 4 ) then
     doubles(n_doubles) = idx(i)
-      n_doubles = n_doubles+1
-    else if ( degree == 2 ) then
     singles(n_singles) = idx(i)
-      n_singles = n_singles+1
+    n_doubles = n_doubles+add_double(degree)
-    endif
+    n_singles = n_singles+add_single(degree)
   enddo
   n_singles = n_singles-1
   n_doubles = n_doubles-1
@@ -1113,15 +1111,14 @@ subroutine get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_
   integer                        :: i
   integer(bit_kind)              :: v
   integer                        :: degree
+  integer                        :: add_single(0:64) = (/ 0, 0, 1, 0, 0, (0, i=1,60) /)
   include 'utils/constants.include.F'
 
   n_singles = 1
   do i=1,size_buffer
     degree = popcnt(xor( spindet, buffer(i) ))
-    if (degree == 2) then
     singles(n_singles) = idx(i)
-      n_singles = n_singles+1
+    n_singles = n_singles+add_single(degree)
-    endif
   enddo
   n_singles = n_singles-1
 
@@ -1145,14 +1142,13 @@ subroutine get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_
   integer                        :: i
   include 'utils/constants.include.F'
   integer                        :: degree
+  integer                        :: add_double(0:64) = (/ 0, 0, 0, 0, 1, (0, i=1,60) /)
 
   n_doubles = 1
   do i=1,size_buffer
     degree = popcnt(xor( spindet, buffer(i) ))
-    if ( degree == 4 ) then
     doubles(n_doubles) = idx(i)
-      n_doubles = n_doubles+1
+    n_doubles = n_doubles+add_double(degree)
-    endif
   enddo
   n_doubles = n_doubles-1
 
@@ -1193,16 +1189,10 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu
       xorvec(k) = xor( spindet(k), buffer(k,i) )
     enddo
 
-    if (xorvec(1) /= 0_8) then
-      degree = popcnt(xorvec(1))
-    else
     degree = 0
-    endif
 
-    do k=2,$N_int
+    do k=1,$N_int
-      if ( (degree <= 4).and.(xorvec(k) /= 0_8) ) then
         degree = degree + popcnt(xorvec(k))
-      endif
     enddo
 
     if ( degree == 4 ) then
@@ -1247,22 +1237,18 @@ subroutine get_all_spin_singles_$N_int(buffer, idx, spindet, size_buffer, single
       xorvec(k) = xor( spindet(k), buffer(k,i) )
     enddo
 
-    if (xorvec(1) /= 0_8) then
-      degree = popcnt(xorvec(1))
-    else
     degree = 0
-    endif
 
-    do k=2,$N_int
+    do k=1,$N_int
-      if ( (degree <= 2).and.(xorvec(k) /= 0_8) ) then
         degree = degree + popcnt(xorvec(k))
-      endif
     enddo
 
-    if ( degree == 2 ) then
+    if ( degree /= 2 ) then
+      cycle
+    endif
+
     singles(n_singles) = idx(i)
     n_singles = n_singles+1
-    endif
 
   enddo
   n_singles = n_singles-1
@@ -1296,22 +1282,18 @@ subroutine get_all_spin_doubles_$N_int(buffer, idx, spindet, size_buffer, double
       xorvec(k) = xor( spindet(k), buffer(k,i) )
     enddo
 
-    if (xorvec(1) /= 0_8) then
-      degree = popcnt(xorvec(1))
-    else
     degree = 0
-    endif
 
-    do k=2,$N_int
+    do k=1,$N_int
-      if ( (degree <= 4).and.(xorvec(k) /= 0_8) ) then
         degree = degree + popcnt(xorvec(k))
-      endif
     enddo
 
-    if ( degree == 4 ) then
+    if ( degree /= 4 ) then
+      cycle
+    endif
+
     doubles(n_doubles) = idx(i)
     n_doubles = n_doubles+1
-    endif
 
   enddo
 

src/tools/print_energy.irp.f

@@ -14,27 +14,5 @@ end
 
 subroutine run
  implicit none
- BEGIN_DOC
+  print *,  psi_energy + nuclear_repulsion
-! E = \sum_{ij} c_i c_j H_{ij}  / \sum{i} c_i^2
- END_DOC
- integer :: i,j
- double precision :: i_H_psi_array(N_states)
- double precision :: E(N_states)
- double precision :: norm(N_states)
-
- E(1:N_states) = nuclear_repulsion
- norm(1:N_states) = 0.d0
- do i=1,N_det
-  call i_H_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, &
-               size(psi_coef,1), N_states, i_H_psi_array)
-  do j=1,N_states
-    norm(j) += psi_coef(i,j)*psi_coef(i,j)
-    E(j) += i_H_psi_array(j) * psi_coef(i,j)
-  enddo
- enddo
-
- print *, 'Energy:'
- do i=1,N_states
-   print *, E(i)/norm(i)
- enddo
 end

src/two_body_rdm/two_e_dm_mo.irp.f

@@ -34,6 +34,7 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num,1)]
      enddo
     enddo
    enddo
+   two_e_dm_mo(:,:,:,:,:) = two_e_dm_mo(:,:,:,:,:) * 2.d0
 
  END_PROVIDER
 

src/utils/c_funcs.c

@@ -1,6 +1,40 @@
 #include <unistd.h>
+#include <stdio.h>
+#include <string.h>
 
 void usleep_c(int s)
 {
   usleep((useconds_t) s);
 }
+
+void sscanf_ssds_c(const char* str, char* s1, char* s2, int* i, char* s3)
+{
+  sscanf(str, "%s %s %d %s", s1, s2, i, s3);
+  s1[strlen(s1)] = ' ';
+  s2[strlen(s2)] = ' ';
+  s3[strlen(s3)] = ' ';
+}
+
+void sscanf_dd_c(const char* str, int* i1, int* i2)
+{
+  sscanf(str, "%d %d", i1, i2);
+}
+
+void sscanf_ddd_c(const char* str, int* i1, int* i2, int* i3)
+{
+  sscanf(str, "%d %d %d", i1, i2, i3);
+}
+
+void sscanf_ss_c(const char* str, char* s1, char* s2)
+{
+  sscanf(str, "%s %s", s1, s2);
+  s1[strlen(s1)] = ' ';
+  s2[strlen(s2)] = ' ';
+}
+
+void sscanf_sd_c(const char* str, char* s1, int* i)
+{
+  sscanf(str, "%s %d", s1, i);
+  s1[strlen(s1)] = ' ';
+}
+

src/utils/c_functions.f90

@@ -8,14 +8,127 @@ module c_functions
      end subroutine usleep_c
   end interface
 
-end module
+  interface
+     integer(c_int) function atoi_c(a) bind (C,name="atoi")
+       use iso_c_binding
+       character(kind=c_char), intent(in) :: a(*)
+     end function atoi_c
+  end interface
 
-subroutine usleep(us)
+  interface
+     subroutine sscanf_ss_c(str,s1, s2) bind (C)
+       use iso_c_binding
+       character(kind=c_char), intent(in ) :: str(*)
+       character(kind=c_char), intent(out) :: s1(*),s2(*)
+     end subroutine sscanf_ss_c
+  end interface
+
+  interface
+     subroutine sscanf_ssds_c(str, s1, s2, i, s3) bind (C)
+       use iso_c_binding
+       character(kind=c_char), intent(in ) :: str(*)
+       character(kind=c_char), intent(out) :: s1(*),s2(*),s3(*)
+       integer(kind=c_int)   , intent(out) :: i
+     end subroutine sscanf_ssds_c
+  end interface
+
+  interface
+     subroutine sscanf_dd_c(str, i1, i2) bind (C)
+       use iso_c_binding
+       character(kind=c_char), intent(in ) :: str(*)
+       integer(kind=c_int)   , intent(out) :: i1, i2
+     end subroutine sscanf_dd_c
+  end interface
+
+  interface
+     subroutine sscanf_ddd_c(str, i1, i2, i3) bind (C)
+       use iso_c_binding
+       character(kind=c_char), intent(in ) :: str(*)
+       integer(kind=c_int)   , intent(out) :: i1, i2, i3
+     end subroutine sscanf_ddd_c
+  end interface
+
+  interface
+     subroutine sscanf_sd_c(str,s1, i) bind (C)
+       use iso_c_binding
+       character(kind=c_char), intent(in ) :: str(*)
+       character(kind=c_char), intent(out) :: s1(*)
+       integer(kind=c_int)   , intent(out) :: i
+     end subroutine sscanf_sd_c
+  end interface
+
+contains
+
+  integer function atoi(a)
+    implicit none
+    character(len=*), intent(in) :: a
+    atoi = atoi_c(trim(a)//c_null_char)
+  end function atoi
+
+end module c_functions
+
+subroutine sscanf_ss(str, s1,s2)
   use c_functions
   use iso_c_binding
   implicit none
+  character(*), intent(in)  :: str
+  character(*), intent(out) :: s1,s2
+  s1 = ' '
+  s2 = ' '
+  call sscanf_ss_c(trim(str)//c_null_char, s1, s2)
+end subroutine sscanf_ss
+
+subroutine sscanf_sd(str, s1,i)
+  use c_functions
+  use iso_c_binding
+  implicit none
+  character(*), intent(in)  :: str
+  character(*), intent(out) :: s1
+  integer, intent(out)      :: i
+  s1 = ' '
+  call sscanf_sd_c(trim(str)//c_null_char, s1, i)
+end subroutine sscanf_sd
+
+subroutine sscanf_ssds(str, s1,s2,i,s3)
+  use c_functions
+  use iso_c_binding
+  implicit none
+  character(*), intent(in)  :: str
+  character(*), intent(out) :: s1,s2,s3
+  integer, intent(out)      :: i
+  s1 = ' '
+  s2 = ' '
+  s3 = ' '
+  call sscanf_ssds_c(trim(str)//c_null_char, s1, s2, i, s3)
+end subroutine sscanf_ssds
+
+subroutine sscanf_dd(str, i1,i2)
+  use c_functions
+  use iso_c_binding
+  implicit none
+  character(*), intent(in)  :: str
+  integer, intent(out)      :: i1, i2
+  call sscanf_dd_c(trim(str)//c_null_char, i1, i2)
+end subroutine sscanf_dd
+
+subroutine sscanf_ddd(str, i1,i2,i3)
+  use c_functions
+  use iso_c_binding
+  implicit none
+  character(*), intent(in)  :: str
+  integer, intent(out)      :: i1, i2, i3
+  call sscanf_ddd_c(trim(str)//c_null_char, i1, i2, i3)
+end subroutine sscanf_ddd
+
+
+subroutine usleep(us)
+  use iso_c_binding
+  use c_functions
+  implicit none
   integer, intent(in) :: us
   integer(c_int) :: u
   u = us
   call usleep_c(u)
-end
+end subroutine usleep
+
+

src/zmq/utils.irp.f

@@ -16,6 +16,7 @@ END_PROVIDER
  BEGIN_PROVIDER [ character*(128), qp_run_address ]
 &BEGIN_PROVIDER [ integer, zmq_port_start ]
   use f77_zmq
+  use c_functions
   implicit none
   BEGIN_DOC
   ! Address of the qp_run socket
@@ -32,14 +33,15 @@ END_PROVIDER
   do i=len(buffer),1,-1
     if ( buffer(i:i) == ':') then
       qp_run_address = trim(buffer(1:i-1))
-      read(buffer(i+1:), *, err=10,end=10) zmq_port_start
+      zmq_port_start = atoi(buffer(i+1:))
       exit
     endif
   enddo
-  return
+
-  10 continue
+  if (zmq_port_start == 0) then
-  print *,  irp_here, ': Error in read'
+    print *,  irp_here, ': zmq_port_start is 0'
     stop -1
+  endif
 END_PROVIDER
 
  BEGIN_PROVIDER [ character*(128), zmq_socket_pull_tcp_address    ]
@@ -84,6 +86,7 @@ end
 
 subroutine switch_qp_run_to_master
   use f77_zmq
+  use c_functions
   implicit none
   BEGIN_DOC
   ! Address of the master qp_run socket
@@ -102,16 +105,17 @@ subroutine switch_qp_run_to_master
   do i=len(buffer),1,-1
     if ( buffer(i:i) == ':') then
       qp_run_address = trim(buffer(1:i-1))
-      read(buffer(i+1:), *, end=10, err=10) zmq_port_start
+      zmq_port_start = atoi(buffer(i+1:))
       exit
     endif
   enddo
   call reset_zmq_addresses
 
   return
-  10 continue
+  if (zmq_port_start == 0) then
-  print *,  irp_here, ': Error in read'
+    print *,  irp_here, ': zmq_port_start is 0'
     stop -1
+  endif
 end
 
 
@@ -650,12 +654,17 @@ integer function connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
   rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 510, 0)
   message = trim(message(1:rc))
   if(message(1:5) == "error") then
-    go to 10
+    connect_to_taskserver = -1
+    return
   end if
-  read(message,*, end=10, err=10) reply, state, worker_id, address
+
+  call sscanf_ssds(message, reply, state, worker_id, address)
+
   if (trim(reply) /= 'connect_reply') then
-    go to 10
+    connect_to_taskserver = -1
+    return
   endif
+
   if (trim(state) /= zmq_state) then
     integer, external :: disconnect_from_taskserver_state
     if (disconnect_from_taskserver_state(zmq_to_qp_run_socket, worker_id, state) == -1) then
@@ -663,13 +672,8 @@ integer function connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
       continue
     endif
     connect_to_taskserver = -1
-    return
   endif
 
-  return
-  10 continue
-!  print *,  irp_here//': '//trim(message)
-  connect_to_taskserver = -1
 end
 
 integer function disconnect_from_taskserver(zmq_to_qp_run_socket, worker_id)
@@ -698,7 +702,7 @@ integer function disconnect_from_taskserver_state(zmq_to_qp_run_socket, worker_i
   character*(512)                :: message, reply
   character*(128)                :: state_tmp
 
-  disconnect_from_taskserver_state = 0
+  disconnect_from_taskserver_state = -1
 
   write(message,*) 'disconnect '//trim(state), worker_id
 
@@ -718,21 +722,15 @@ integer function disconnect_from_taskserver_state(zmq_to_qp_run_socket, worker_i
   rc = min(510,rc)
   message = trim(message(1:rc))
 
-  read(message,*, end=10, err=10) reply, state_tmp
+  call sscanf_ss(message, reply, state_tmp)
-  if ((trim(reply) == 'disconnect_reply').and.(trim(state_tmp) == trim(state))) then
+
-    return
+  if (trim(state_tmp) /= trim(state)) then
-  endif
-  if (trim(message) == 'error Wrong state') then
-    disconnect_from_taskserver_state = -1
-    return
-  else if (trim(message) == 'error No job is running') then
-    disconnect_from_taskserver_state = -1
     return
   endif
 
-  return
+  if ((trim(reply) == 'disconnect_reply')) then
-  10 continue
+    disconnect_from_taskserver_state = 0
-  disconnect_from_taskserver_state = -1
+  endif
 end
 
 integer function add_task_to_taskserver(zmq_to_qp_run_socket,task)
@@ -898,7 +896,7 @@ integer function get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id
 
   character*(1024)                :: message
   character*(64)                 :: reply
-  integer                        :: rc, sze
+  integer                        :: rc, sze, i
 
   get_task_from_taskserver = 0
 
@@ -912,12 +910,15 @@ integer function get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id
   endif
 
   task_id = 0
-  message = repeat(' ',1024)
+  message = ' '
   rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 1024, 0)
-  rc = min(64,rc)
+  i = 1
-  read(message(1:rc),*, end=10, err=10) reply
+  do while (message(i:i) /= ' ')
-  if (trim(reply) == 'get_task_reply') then
+    i = i+1
-    read(message(1:rc),*, end=10, err=10) reply, task_id
+  enddo
+  reply = message(1:i-1)
+  if (reply == 'get_task_reply') then
+    call sscanf_sd(message, reply, task_id)
     rc = 15
     do while (rc < 1024 .and. message(rc:rc) == ' ')
       rc += 1
@@ -937,15 +938,12 @@ integer function get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id
     get_task_from_taskserver = -1
     return
   endif
-  return
-
-  10 continue
-  get_task_from_taskserver = -1
 
 end
 
 
 integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id,task,n_tasks)
+  use c_functions
   use f77_zmq
   implicit none
   BEGIN_DOC
@@ -1000,7 +998,7 @@ integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_i
       return
     endif
     rc = min(1024,rc)
-    read(message(1:rc),*, end=10, err=10) task_id(i)
+    task_id(i) = atoi(message(1:rc))
     if (task_id(i) == 0) then
       task(i) = 'terminate'
       n_tasks = i