diff --git a/ocaml/Makefile b/ocaml/Makefile
index b666187f..a49304e6 100644
--- a/ocaml/Makefile
+++ b/ocaml/Makefile
@@ -12,7 +12,8 @@ endif
 LIBS=
 PKGS=
 OCAMLCFLAGS="-g"
-OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) 
+OCAMLOPTFLAGS="opt -O3 -nodynlink -remove-unused-arguments -rounds 16 -inline 100 -inline-max-unroll 100"
+OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) -ocamlopt $(OCAMLOPTFLAGS)
 MLLFILES=$(wildcard *.mll) 
 MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml Input_auto_generated.ml qp_edit.ml 
 MLIFILES=$(wildcard *.mli) git
diff --git a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
index 8d36fd55..04718218 100644
--- a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
+++ b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
@@ -396,7 +396,11 @@ end function
 
 BEGIN_PROVIDER [ integer, comb_teeth ]
   implicit none
-  comb_teeth = 100
+  BEGIN_DOC
+! Number of teeth in the comb
+  END_DOC
+  comb_teeth = min(1+N_det/10,100)
+
 END_PROVIDER
 
 
@@ -555,7 +559,7 @@ END_PROVIDER
   comb_step = 1d0/dfloat(comb_teeth)
   first_det_of_comb = 1
   do i=1,N_det_generators
-    if(pt2_weight(i)/norm_left < .5d0*comb_step) then
+    if(pt2_weight(i)/norm_left < .25d0*comb_step) then
       first_det_of_comb = i
       exit
     end if
diff --git a/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f b/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f
index 613f00e5..4d5bb344 100644
--- a/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f
+++ b/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f
@@ -82,8 +82,8 @@ subroutine run_pt2_slave(thread,iproc,energy)
     endif
     call push_pt2_results(zmq_socket_push, i_generator, pt2, task_id, n_tasks)
 
-    ! Try to adjust n_tasks around 10 seconds per job
-    n_tasks = int(10.d0 * dble(n_tasks) / (time1 - time0 + 1.d-9))+1
+    ! Try to adjust n_tasks around 1 second per job
+    n_tasks = min(n_tasks,int( dble(n_tasks) / (time1 - time0 + 1.d-9)))+1
   end do
 
   integer, external :: disconnect_from_taskserver
diff --git a/plugins/Full_CI_ZMQ/run_selection_slave.irp.f b/plugins/Full_CI_ZMQ/run_selection_slave.irp.f
index 43cd340a..6de5b19d 100644
--- a/plugins/Full_CI_ZMQ/run_selection_slave.irp.f
+++ b/plugins/Full_CI_ZMQ/run_selection_slave.irp.f
@@ -97,8 +97,8 @@ subroutine run_selection_slave_new(thread,iproc,energy)
     pt2(:,:) = 0d0
     buf%cur = 0
 
-    ! Try to adjust n_tasks around 10 seconds per job
-    n_tasks = int(10.d0 * dble(n_tasks) / (time1 - time0 + 1.d-9))+1
+    ! Try to adjust n_tasks around 1 second per job
+    n_tasks = min(n_tasks,int( dble(n_tasks) / (time1 - time0 + 1.d-9)))+1
   end do
 
   integer, external :: disconnect_from_taskserver
diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f
index 3da1bca4..ba02dbf9 100644
--- a/plugins/Full_CI_ZMQ/selection.irp.f
+++ b/plugins/Full_CI_ZMQ/selection.irp.f
@@ -5,7 +5,7 @@ BEGIN_PROVIDER [ integer, fragment_count ]
   BEGIN_DOC
   ! Number of fragments for the deterministic part
   END_DOC
-  fragment_count = (elec_alpha_num-n_core_orb)*mo_tot_num
+  fragment_count = (elec_alpha_num-n_core_orb)**2
 END_PROVIDER
 
 
@@ -690,8 +690,8 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
     
     if (interesting(i) == i_gen) then
         if(sp == 3) then
-          do j=1,mo_tot_num
-            do k=1,mo_tot_num
+          do k=1,mo_tot_num
+            do j=1,mo_tot_num
               banned(j,k,2) = banned(k,j,1)
             enddo
           enddo
@@ -909,14 +909,22 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
     if(.not. bannedOrb(puti, mi)) then
       tmp_row = 0d0
       do putj=1, hfix-1
-        if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
-        hij = (mo_bielec_integral(putj, hfix, p1, p2)-mo_bielec_integral(putj,hfix,p2,p1)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
-        tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
+        if(lbanned(putj, ma)) cycle
+        if(banned(putj, puti,bant)) cycle
+        hij = mo_bielec_integral(putj, hfix, p1, p2)-mo_bielec_integral(putj,hfix,p2,p1)
+        if (hij /= 0.d0) then
+          hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
+          tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
+        endif
       end do
       do putj=hfix+1, mo_tot_num
-        if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
-        hij = (mo_bielec_integral(hfix,putj,p1, p2)-mo_bielec_integral(hfix,putj,p2,p1)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
-        tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
+        if(lbanned(putj, ma)) cycle
+        if(banned(putj, puti,bant)) cycle
+        hij = mo_bielec_integral(putj,hfix,p2, p1)-mo_bielec_integral(putj,hfix,p1,p2)
+        if (hij /= 0.d0) then
+          hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
+          tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
+        endif
       end do
 
       if(ma == 1) then           
@@ -935,14 +943,20 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
       !p1 fixed
       putj = p1
       if(.not. banned(putj,puti,bant)) then
-        hij = mo_bielec_integral(puti,hfix,pfix,p2) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix)
-        tmp_row(:,puti) += hij * coefs(:)
+        hij = mo_bielec_integral(puti,hfix,pfix,p2)
+        if (hij /= 0.d0) then
+          hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix)
+          tmp_row(:,puti) += hij * coefs(:)
+        endif
       end if
       
       putj = p2
       if(.not. banned(putj,puti,bant)) then
-        hij = mo_bielec_integral(puti,hfix,pfix,p1) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix)
-        tmp_row2(:,puti) += hij * coefs(:)
+        hij = mo_bielec_integral(puti,hfix,pfix,p1)
+        if (hij /= 0.d0) then
+          hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix)
+          tmp_row2(:,puti) += hij * coefs(:)
+        endif
       end if
     end do
     
@@ -953,7 +967,9 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
       mat(:,p1,:) += tmp_row(:,:)
       mat(:,p2,:) += tmp_row2(:,:)
     end if
-  else
+
+  else  ! sp /= 3
+
     if(p(0,ma) == 3) then
       do i=1,3
         hfix = h(1,ma)
@@ -962,14 +978,24 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         p2 = p(turn3(2,i), ma)
         tmp_row = 0d0
         do putj=1,hfix-1
-          if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
-          hij = (mo_bielec_integral(putj, hfix, p1, p2)-mo_bielec_integral(putj,hfix,p2,p1)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
-          tmp_row(:,putj) += hij * coefs(:)
+          if(lbanned(putj,ma)) cycle
+!          if(banned(puti,putj,1)) cycle
+          if(banned(putj,puti,1)) cycle
+          hij = mo_bielec_integral(putj, hfix, p1, p2)-mo_bielec_integral(putj,hfix,p2,p1)
+          if (hij /= 0.d0) then
+            hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
+            tmp_row(:,putj) += hij * coefs(:)
+          endif
         end do
         do putj=hfix+1,mo_tot_num
-          if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
-          hij = (mo_bielec_integral(putj, hfix, p2, p1)-mo_bielec_integral(putj,hfix,p1,p2)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
-          tmp_row(:,putj) += hij * coefs(:)
+          if(lbanned(putj,ma)) cycle
+!          if(banned(puti,putj,1)) cycle
+          if(banned(putj,puti,1)) cycle
+          hij = mo_bielec_integral(putj, hfix, p2, p1)-mo_bielec_integral(putj,hfix,p1,p2)
+          if (hij /= 0.d0) then
+            hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
+            tmp_row(:,putj) += hij * coefs(:)
+          endif
         end do
 
         mat(:, :puti-1, puti) += tmp_row(:,:puti-1)
@@ -986,14 +1012,20 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
         if(lbanned(puti,ma)) cycle
         putj = p2
         if(.not. banned(puti,putj,1)) then
-          hij = mo_bielec_integral(puti,hfix, p1,pfix) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1)
-          tmp_row(:,puti) += hij * coefs(:)
+          hij = mo_bielec_integral(puti,hfix, p1,pfix)
+          if (hij /= 0.d0) then
+            hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1)
+            tmp_row(:,puti) += hij * coefs(:)
+          endif
         end if
         
         putj = p1
         if(.not. banned(puti,putj,1)) then
-          hij = mo_bielec_integral(puti, hfix, p2, pfix) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2)
-          tmp_row2(:,puti) += hij * coefs(:)
+          hij = mo_bielec_integral(puti, hfix, p2, pfix) 
+          if (hij /= 0.d0) then
+            hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2)
+            tmp_row2(:,puti) += hij * coefs(:)
+          endif
         end if
       end do
       mat(:,:p2-1,p2) += tmp_row(:,:p2-1)
@@ -1016,10 +1048,11 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
     do i1=1,p(0,s1)
       ib = 1
       if(s1 == s2) ib = i1+1
+      p1 = p(i1,s1)
+      if(bannedOrb(p1, s1)) cycle
       do i2=ib,p(0,s2)
-        p1 = p(i1,s1)
         p2 = p(i2,s2)
-        if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
+        if(bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
         call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
         call i_h_j(gen, det, N_int, hij)
         mat(:, p1, p2) += coefs(:) * hij
diff --git a/plugins/Full_CI_ZMQ/zmq_selection.irp.f b/plugins/Full_CI_ZMQ/zmq_selection.irp.f
index 60d0dbe1..55cdb718 100644
--- a/plugins/Full_CI_ZMQ/zmq_selection.irp.f
+++ b/plugins/Full_CI_ZMQ/zmq_selection.irp.f
@@ -60,7 +60,7 @@ subroutine ZMQ_selection(N_in, pt2)
   task = ' ' 
 
   do i= 1, N_det_generators
-    if (i>ishft(N_det_generators,-7)) then
+    if (i>ishft(N_det_generators,-2)) then
       write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') 0, i, N
       ipos += 30
       if (ipos > 63970) then
diff --git a/src/Integrals_Bielec/qp_ao_ints.irp.f b/src/Integrals_Bielec/qp_ao_ints.irp.f
index 93f62a7d..f499c848 100644
--- a/src/Integrals_Bielec/qp_ao_ints.irp.f
+++ b/src/Integrals_Bielec/qp_ao_ints.irp.f
@@ -1,10 +1,16 @@
 program qp_ao_ints
   use omp_lib
   implicit none
+  IRP_IF MPI
+    include 'mpif.h'
+  IRP_ENDIF
+  integer :: ierr
+
   BEGIN_DOC
 ! Increments a running calculation to compute AO integrals
   END_DOC
   integer :: i
+  PROVIDE zmq_context mpi_master zmq_state zmq_context
 
   call switch_qp_run_to_master
 
@@ -17,15 +23,27 @@ program qp_ao_ints
   double precision :: integral, ao_bielec_integral
   integral = ao_bielec_integral(1,1,1,1)
 
-  character*(64) :: state
-  call wait_for_state(zmq_state,state)
-  do while (state /= 'Stopped')
+  do
+    call wait_for_state('ao_integrals',zmq_state)
+    if (zmq_state(1:7) == 'Stopped') then
+      exit
+    endif
+
     !$OMP PARALLEL DEFAULT(PRIVATE) PRIVATE(i)
     i = omp_get_thread_num()
     call ao_bielec_integrals_in_map_slave_tcp(i)
     !$OMP END PARALLEL
-    call wait_for_state(zmq_state,state)
+    IRP_IF MPI
+      call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+      if (ierr /= MPI_SUCCESS) then
+        print *,  irp_here, 'error in barrier'
+      endif
+    IRP_ENDIF
+
   enddo
+  IRP_IF MPI
+    call MPI_finalize(i)
+  IRP_ENDIF
 
   print *,  'Done'
 end
diff --git a/src/ZMQ/utils.irp.f b/src/ZMQ/utils.irp.f
index f8f7c71b..c263c884 100644
--- a/src/ZMQ/utils.irp.f
+++ b/src/ZMQ/utils.irp.f
@@ -337,7 +337,7 @@ IRP_ENDIF
     stop 'Unable to set ZMQ_LINGER on push socket'
   endif
   
-  rc = f77_zmq_setsockopt(new_zmq_push_socket,ZMQ_SNDHWM,5,4)
+  rc = f77_zmq_setsockopt(new_zmq_push_socket,ZMQ_SNDHWM,1,4)
   if (rc /= 0) then
     stop 'Unable to set ZMQ_SNDHWM on push socket'
   endif