diff --git a/plugins/QMC/EZFIO.cfg b/plugins/QMC/EZFIO.cfg
new file mode 100644
index 00000000..8960a0fa
--- /dev/null
+++ b/plugins/QMC/EZFIO.cfg
@@ -0,0 +1,6 @@
+[ci_threshold]
+type: Threshold
+doc: Threshold on the CI coefficients as in QMCChem
+interface: ezfio,provider,ocaml
+default: 1.e-8
+
diff --git a/plugins/QMC/NEEDED_CHILDREN_MODULES b/plugins/QMC/NEEDED_CHILDREN_MODULES
index 34de8ddb..9a2f60c0 100644
--- a/plugins/QMC/NEEDED_CHILDREN_MODULES
+++ b/plugins/QMC/NEEDED_CHILDREN_MODULES
@@ -1 +1 @@
-Determinants Davidson
+Determinants Davidson Full_CI_ZMQ
diff --git a/plugins/QMC/save_for_qmcchem.irp.f b/plugins/QMC/save_for_qmcchem.irp.f
index a281a184..771bf618 100644
--- a/plugins/QMC/save_for_qmcchem.irp.f
+++ b/plugins/QMC/save_for_qmcchem.irp.f
@@ -24,13 +24,13 @@ program save_for_qmc
   )
   iunit = 13
   open(unit=iunit,file=trim(ezfio_filename)//'/simulation/e_ref',action='write')
-  call ezfio_has_full_ci_energy_pt2(exists)
+  call ezfio_has_full_ci_zmq_energy_pt2(exists)
   if (exists) then
-    call ezfio_get_full_ci_energy_pt2(e_ref)
+    call ezfio_get_full_ci_zmq_energy_pt2(e_ref)
   else
-    call ezfio_has_full_ci_energy(exists)
+    call ezfio_has_full_ci_zmq_energy(exists)
     if (exists) then
-      call ezfio_get_full_ci_energy(e_ref)
+      call ezfio_get_full_ci_zmq_energy(e_ref)
     else
       call ezfio_has_hartree_fock_energy(exists)
       if (exists) then
diff --git a/plugins/QMC/save_for_qmcpack.irp.f b/plugins/QMC/save_for_qmcpack.irp.f
index 186ca616..51ed866d 100644
--- a/plugins/QMC/save_for_qmcpack.irp.f
+++ b/plugins/QMC/save_for_qmcpack.irp.f
@@ -21,6 +21,6 @@ program qmcpack
   enddo
   call save_mos
   call system('rm '//trim(ezfio_filename)//'/mo_basis/ao_md5')
-  call system('$QP_ROOT/src/qmcpack/qp_convert_qmcpack_to_ezfio.py '//trim(ezfio_filename))
+  call system('$QP_ROOT/src/QMC/qp_convert_qmcpack_to_ezfio.py '//trim(ezfio_filename))
 
 end
diff --git a/plugins/QMC/truncate_wf_spin.irp.f b/plugins/QMC/truncate_wf_spin.irp.f
index b4d6e500..5a5fe125 100644
--- a/plugins/QMC/truncate_wf_spin.irp.f
+++ b/plugins/QMC/truncate_wf_spin.irp.f
@@ -1,4 +1,10 @@
-program e_curve
+program truncate
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call run
+end
+
+subroutine run
   use bitmasks
   implicit none
   integer :: i,j,k, kk, nab, m, l
@@ -6,24 +12,17 @@ program e_curve
   integer, allocatable :: iorder(:)
   double precision , allocatable :: norm_sort(:)
   double precision :: e_0(N_states)
-  if (.not.read_wf) then
-    stop 'Please set read_wf to true'
-  endif
 
   PROVIDE mo_bielec_integrals_in_map H_apply_buffer_allocated
 
   nab = n_det_alpha_unique+n_det_beta_unique
   allocate ( norm_sort(0:nab), iorder(0:nab) )
 
-  double precision :: thresh
   integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
   double precision, allocatable  :: u_t(:,:), v_t(:,:), s_t(:,:)
   double precision, allocatable  :: u_0(:,:), v_0(:,:)
   allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det))
-  allocate(u_0(N_states,N_det),v_0(N_states,N_det))
-
-  print *,  'Threshold?'
-  read(*,*) thresh
+  allocate(u_0(N_det,N_states),v_0(N_det,N_states))
 
   norm_sort(0) = 0.d0
   iorder(0) = 0
@@ -45,19 +44,23 @@ program e_curve
   do j=0,nab
     i = iorder(j)
     if (i<0) then
+      !$OMP PARALLEL DO PRIVATE(k)
       do k=1,n_det
         if (psi_bilinear_matrix_columns(k) == -i) then
           psi_bilinear_matrix_values(k,1) = 0.d0
         endif
       enddo
+      !$OMP END PARALLEL DO
     else
+      !$OMP PARALLEL DO PRIVATE(k)
       do k=1,n_det
         if (psi_bilinear_matrix_rows(k) ==  i) then
           psi_bilinear_matrix_values(k,1) = 0.d0
         endif
       enddo
+      !$OMP END PARALLEL DO
     endif
-    if (thresh > norm_sort(j)) then
+    if (ci_threshold > norm_sort(j)) then
       cycle
     endif
 
@@ -70,7 +73,9 @@ program e_curve
         u_t,                                                           &
         size(u_t, 1),                                                  &
         N_det, N_states)
+    print *, 'Computing H|Psi> ...'
     call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
+    print *, 'Done'
     call dtranspose(                                                   &
         v_t,                                                           &
         size(v_t, 1),                                                  &
@@ -96,7 +101,7 @@ program e_curve
     print *,  'Energy', E
     exit
   enddo
-  call wf_of_psi_bilinear_matrix()
+  call wf_of_psi_bilinear_matrix(.True.)
   call save_wavefunction
 
   deallocate (iorder, norm_sort)
diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f
index 51572462..e3f5c0b1 100644
--- a/src/Determinants/slater_rules.irp.f
+++ b/src/Determinants/slater_rules.irp.f
@@ -43,10 +43,12 @@ subroutine get_excitation_degree(key1,key2,degree,Nint)
       degree = sum(popcnt(xorvec(1:8)))
 
     case default
-      do l=1,ishft(Nint,1)
+      integer :: lmax
+      lmax = ishft(Nint,1)
+      do l=1,lmax
         xorvec(l) = xor( key1(l), key2(l))
       enddo
-      degree = sum(popcnt(xorvec(1:l)))
+      degree = sum(popcnt(xorvec(1:lmax)))
   
   end select
 
@@ -245,12 +247,16 @@ subroutine get_double_excitation(det1,det2,exc,phase,Nint)
         if (j==k) then
           nperm = nperm + popcnt(iand(det1(j,ispin),                 &
               iand( ibset(0_bit_kind,m-1)-1_bit_kind,                &
-              ibclr(-1_bit_kind,n)+1_bit_kind ) ))
+                    ibclr(-1_bit_kind,n)+1_bit_kind ) ))
+! TODO        iand( not(ishft(1_bit_kind,n+1))+1_bit_kind, &
+!                   ishft(1_bit_kind,m)-1_bit_kind)))
         else
           nperm = nperm + popcnt(iand(det1(k,ispin),                 &
               ibset(0_bit_kind,m-1)-1_bit_kind)) 
+! TODO        ishft(1_bit_kind,m)-1_bit_kind)) 
           if (n < bit_kind_size) then
               nperm = nperm + popcnt(iand(det1(j,ispin), ibclr(-1_bit_kind,n) +1_bit_kind))
+! TODO        ishft(1_bit_kind,m)-1_bit_kind)) 
           endif
           do i=j+1,k-1
             nperm = nperm + popcnt(det1(i,ispin))
@@ -365,8 +371,12 @@ subroutine get_mono_excitation(det1,det2,exc,phase,Nint)
       if (j==k) then
         nperm = popcnt(iand(det1(j,ispin),                           &
             iand(ibset(0_bit_kind,m-1)-1_bit_kind,ibclr(-1_bit_kind,n)+1_bit_kind)))
+!TODO       iand( not(ishft(1_bit_kind,n+1))+1_bit_kind, &
+!                 ishft(1_bit_kind,m)-1_bit_kind)))
       else
         nperm = nperm + popcnt(iand(det1(k,ispin),ibset(0_bit_kind,m-1)-1_bit_kind))
+!TODO    nperm = popcnt(iand(det1(k,ispin), ishft(1_bit_kind,m)-1_bit_kind)) + &
+!                popcnt(iand(det1(j,ispin), not(ishft(1_bit_kind,n+1))+1_bit_kind))
         if (n < bit_kind_size) then
             nperm = nperm + popcnt(iand(det1(j,ispin),ibclr(-1_bit_kind,n)+1_bit_kind))
         endif
@@ -1495,14 +1505,16 @@ subroutine get_excitation_degree_vector_double_alpha_beta(key1,key2,degree,Nint,
     !DIR$ LOOP COUNT (1000)
     do i=1,sze
       d = 0
+      degree_alpha = 0
+      degree_beta  = 0
       !DIR$ LOOP COUNT MIN(4)
       do m=1,Nint
         d = d + popcnt(xor( key1(m,1,i), key2(m,1)))                 &
               + popcnt(xor( key1(m,2,i), key2(m,2)))
         key_tmp(m,1) = xor(key1(m,1,i),key2(m,1))
         key_tmp(m,2) = xor(key1(m,2,i),key2(m,2))
-        degree_alpha = popcnt(key_tmp(m,1)) 
-        degree_beta  = popcnt(key_tmp(m,2)) 
+        degree_alpha += popcnt(key_tmp(m,1)) 
+        degree_beta  += popcnt(key_tmp(m,2)) 
       enddo
       if(degree_alpha .gt.3 .or. degree_beta .gt.3 )cycle !! no double excitations of same spin
         degree(l) = ishft(d,-1)
@@ -1653,12 +1665,13 @@ subroutine get_excitation_degree_vector_mono_or_exchange_verbose(key1,key2,degre
     do i=1,sze
       d = 0
       exchange_1 = 0
+      exchange_2 = 0
       !DIR$ LOOP COUNT MIN(4)
       do m=1,Nint
         d = d + popcnt(xor( key1(m,1,i), key2(m,1)))                 &
               + popcnt(xor( key1(m,2,i), key2(m,2)))
-        exchange_1 = popcnt(xor(iand(key1(m,1,i),key1(m,2,i)),iand(key2(m,1),key2(m,2))))  
-        exchange_2 = popcnt(iand(xor(key1(m,1,i),key2(m,1)),xor(key1(m,2,i),key2(m,2))))   
+        exchange_1 += popcnt(xor(iand(key1(m,1,i),key1(m,2,i)),iand(key2(m,1),key2(m,2))))  
+        exchange_2 += popcnt(iand(xor(key1(m,1,i),key2(m,1)),xor(key1(m,2,i),key2(m,2))))   
       enddo
       if (d > 4)cycle
       if (d ==4)then  
@@ -2217,7 +2230,7 @@ subroutine get_excitation_degree_spin(key1,key2,degree,Nint)
       degree = sum(popcnt(xorvec(1:4)))
 
     case default
-      do l=1,N_int
+      do l=1,Nint
         xorvec(l) = xor( key1(l), key2(l))
       enddo
       degree = sum(popcnt(xorvec(1:Nint)))
diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f
index a1ad04fe..b6ca1cba 100644
--- a/src/Determinants/spindeterminants.irp.f
+++ b/src/Determinants/spindeterminants.irp.f
@@ -1217,7 +1217,6 @@ subroutine wf_of_psi_bilinear_matrix(truncate)
   integer                        :: idx
   integer, external              :: get_index_in_psi_det_sorted_bit
   double precision               :: norm(N_states)
-  PROVIDE psi_bilinear_matrix
 
   do k=1,N_det
    i = psi_bilinear_matrix_rows(k)
diff --git a/src/Integrals_Bielec/four_idx_transform.irp.f b/src/Integrals_Bielec/four_idx_transform.irp.f
new file mode 100644
index 00000000..a98ce2aa
--- /dev/null
+++ b/src/Integrals_Bielec/four_idx_transform.irp.f
@@ -0,0 +1,12 @@
+program four_idx
+  implicit none
+  BEGIN_DOC
+! 4-index transformation from AO to MO integrals
+  END_DOC
+  
+  disk_access_mo_integrals = 'Write'
+  SOFT_TOUCH disk_access_mo_integrals
+  if (.true.) then
+    PROVIDE mo_bielec_integrals_in_map
+  endif
+end