diff --git a/plugins/Bk/EZFIO.cfg b/plugins/Bk/EZFIO.cfg
new file mode 100644
index 00000000..52d41568
--- /dev/null
+++ b/plugins/Bk/EZFIO.cfg
@@ -0,0 +1,17 @@
+[energy]
+type: double precision
+doc: Calculated energy
+interface: ezfio
+
+[thresh_dressed_ci]
+type: Threshold
+doc: Threshold on the convergence of the dressed CI energy
+interface: ezfio,provider,ocaml
+default: 1.e-5
+
+[n_it_max_dressed_ci]
+type: Strictly_positive_int
+doc: Maximum number of dressed CI iterations
+interface: ezfio,provider,ocaml
+default: 10
+
diff --git a/plugins/Bk/NEEDED_CHILDREN_MODULES b/plugins/Bk/NEEDED_CHILDREN_MODULES
new file mode 100644
index 00000000..6bcca9aa
--- /dev/null
+++ b/plugins/Bk/NEEDED_CHILDREN_MODULES
@@ -0,0 +1,2 @@
+Bitmask dress_zmq DavidsonDressed
+
diff --git a/plugins/Bk/README.rst b/plugins/Bk/README.rst
new file mode 100644
index 00000000..7b379f8e
--- /dev/null
+++ b/plugins/Bk/README.rst
@@ -0,0 +1,12 @@
+==
+Bk
+==
+
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
diff --git a/plugins/Bk/bk.irp.f b/plugins/Bk/bk.irp.f
new file mode 100644
index 00000000..254320f7
--- /dev/null
+++ b/plugins/Bk/bk.irp.f
@@ -0,0 +1,26 @@
+program bk
+  implicit none
+  BEGIN_DOC
+! Shifted-Bk method
+  END_DOC
+  read_wf = .True.
+  state_following = .True.
+  TOUCH read_wf state_following
+  call run()
+end
+
+subroutine run
+  implicit none
+  call diagonalize_ci_dressed
+  integer :: istate
+  print *,  'Bk Energy'
+  print *,  '---------'
+  print *,  ''
+  do istate = 1,N_states
+   print *,  istate, CI_energy_dressed(istate)
+  enddo
+!  call save_wavefunction
+  call ezfio_set_bk_energy(ci_energy_dressed(1))
+end
+
+
diff --git a/plugins/Bk/dressing.irp.f b/plugins/Bk/dressing.irp.f
new file mode 100644
index 00000000..ffd308fa
--- /dev/null
+++ b/plugins/Bk/dressing.irp.f
@@ -0,0 +1,46 @@
+subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, minilist, det_minilist, n_minilist, alpha, iproc)
+ use bitmasks
+ implicit none
+  BEGIN_DOC
+  !delta_ij_loc(:,:,1) : dressing column for H
+  !delta_ij_loc(:,:,2) : dressing column for S2
+  !minilist : indices of determinants connected to alpha ( in psi_det_sorted )
+  !n_minilist : size of minilist
+  !alpha : alpha determinant
+  END_DOC
+  integer, intent(in)             :: Nint, Ndet, Nstates, n_minilist, iproc
+  integer(bit_kind), intent(in)   :: alpha(Nint,2), det_minilist(Nint, 2, n_minilist)
+  integer,intent(in)              :: minilist(n_minilist)
+  double precision, intent(inout) :: delta_ij_loc(Nstates,Ndet,2)
+
+  integer :: j, j_mini, i_state
+  double precision :: c_alpha(N_states), h_alpha_alpha, hdress, sdress
+  double precision :: i_h_alpha, i_s_alpha, alpha_h_psi(N_states)
+
+  double precision, external :: diag_H_mat_elem
+
+  h_alpha_alpha = diag_h_mat_elem(alpha,N_int)
+  call i_H_psi_minilist(alpha,det_minilist,minilist,n_minilist,psi_coef,N_int,n_minilist,size(psi_coef,1),N_states,alpha_h_psi)
+
+  do i_state=1,N_states
+    c_alpha(i_state) = alpha_h_psi(i_state) / &
+      (dress_e0_denominator(i_state) - h_alpha_alpha) 
+  enddo
+
+  do j_mini=1,n_minilist
+    j = minilist(j_mini)
+    call i_H_j (det_minilist(1,1,j_mini),alpha,N_int,i_h_alpha)
+    call get_s2(det_minilist(1,1,j_mini),alpha,N_int,i_s_alpha)
+    do i_state=1,N_states
+      hdress = c_alpha(i_state) * i_h_alpha
+      sdress = c_alpha(i_state) * i_s_alpha
+      delta_ij_loc(i_state,j,1) = delta_ij_loc(i_state,j,1) + hdress
+      delta_ij_loc(i_state,j,2) = delta_ij_loc(i_state,j,2) + sdress
+    enddo
+  enddo
+
+
+end subroutine
+
+
+
diff --git a/plugins/DavidsonDressed/diagonalize_CI.irp.f b/plugins/DavidsonDressed/diagonalize_CI.irp.f
index 7b12bc1c..3d1c1118 100644
--- a/plugins/DavidsonDressed/diagonalize_CI.irp.f
+++ b/plugins/DavidsonDressed/diagonalize_CI.irp.f
@@ -90,81 +90,31 @@ END_PROVIDER
     
     allocate (eigenvectors(size(H_matrix_dressed,1),N_det))
     allocate (eigenvalues(N_det))
-    call lapack_diag(eigenvalues,eigenvectors,                         &
-        H_matrix_dressed,size(H_matrix_dressed,1),N_det)
-    CI_electronic_energy_dressed(:) = 0.d0
-    if (s2_eig) then
-      i_state = 0
-      allocate (s2_eigvalues(N_det))
-      allocate(index_good_state_array(N_det),good_state_array(N_det))
-      good_state_array = .False.
-      call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,   &
-          N_det,size(eigenvectors,1))
-      do j=1,N_det
-        ! Select at least n_states states with S^2 values closed to "expected_s2"
-        if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
-          i_state +=1
-          index_good_state_array(i_state) = j
-          good_state_array(j) = .True.
-        endif
-        if(i_state.eq.N_states) then
-          exit
-        endif
+
+    do j=1,min(N_states,N_det)
+      do i=1,N_det
+        eigenvectors(i,j) = psi_coef(i,j)
       enddo
-      if(i_state .ne.0)then
-        ! Fill the first "i_state" states that have a correct S^2 value
-        do j = 1, i_state
-          do i=1,N_det
-            CI_eigenvectors_dressed(i,j) = eigenvectors(i,index_good_state_array(j))
-          enddo
-          CI_electronic_energy_dressed(j) = eigenvalues(index_good_state_array(j))
-          CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(j))
-        enddo
-        i_other_state = 0
-        do j = 1, N_det
-          if(good_state_array(j))cycle
-          i_other_state +=1
-          if(i_state+i_other_state.gt.n_states_diag)then
-            exit
-          endif
-          do i=1,N_det
-            CI_eigenvectors_dressed(i,i_state+i_other_state) = eigenvectors(i,j)
-          enddo
-          CI_electronic_energy_dressed(i_state+i_other_state) = eigenvalues(j)
-          CI_eigenvectors_s2_dressed(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state)
-        enddo
-        
-      else
-        print*,''
-        print*,'!!!!!!!!   WARNING  !!!!!!!!!'
-        print*,'  Within the ',N_det,'determinants selected'
-        print*,'  and the ',N_states_diag,'states requested'
-        print*,'  We did not find any state with S^2 values close to ',expected_s2
-        print*,'  We will then set the first N_states eigenvectors of the H matrix'
-        print*,'  as the CI_eigenvectors_dressed'
-        print*,'  You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
-        print*,''
-        do j=1,min(N_states_diag,N_det)
-          do i=1,N_det
-            CI_eigenvectors_dressed(i,j) = eigenvectors(i,j)
-          enddo
-          CI_electronic_energy_dressed(j) = eigenvalues(j)
-          CI_eigenvectors_s2_dressed(j) = s2_eigvalues(j)
-        enddo
-      endif
-      deallocate(index_good_state_array,good_state_array)
-      deallocate(s2_eigvalues)
-    else
-      call u_0_S2_u_0(CI_eigenvectors_s2_dressed,eigenvectors,N_det,psi_det,N_int,&
-          min(N_det,N_states_diag),size(eigenvectors,1))
-      ! Select the "N_states_diag" states of lowest energy
-      do j=1,min(N_det,N_states_diag)
+    enddo
+    do mrcc_state=1,N_states
+      do j=mrcc_state,min(N_states,N_det)
         do i=1,N_det
-          CI_eigenvectors_dressed(i,j) = eigenvectors(i,j)
+          eigenvectors(i,j) = psi_coef(i,j)
         enddo
-        CI_electronic_energy_dressed(j) = eigenvalues(j)
       enddo
-    endif
+
+    call lapack_diag(eigenvalues,eigenvectors,                         &
+        H_matrix_dressed(1,1,mrcc_state),size(H_matrix_dressed,1),N_det)
+      CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state)
+      CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state)
+    enddo
+    do k=N_states+1,N_states_diag
+      CI_eigenvectors_dressed(1:N_det,k) = eigenvectors(1:N_det,k)
+      CI_electronic_energy_dressed(k) = eigenvalues(k)
+    enddo
+    call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
+        N_states_diag,size(CI_eigenvectors_dressed,1))
+
     deallocate(eigenvectors,eigenvalues)
   endif
 
@@ -192,19 +142,19 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
  BEGIN_DOC
  ! Dressed H with Delta_ij
  END_DOC
- integer                        :: i, j,istate,ii,jj
- do istate = 1,N_states
+ integer                        :: i, j, ii,jj, dressing_state
+ do dressing_state = 1,N_states
    do j=1,N_det
      do i=1,N_det
-       h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) 
+       h_matrix_dressed(i,j,dressing_state) = h_matrix_all_dets(i,j) 
      enddo
    enddo
-   i = dressed_column_idx(istate)
+   i = dressed_column_idx(dressing_state)
    do j = 1, N_det
-     h_matrix_dressed(i,j,istate) += dressing_column_h(j,istate)
-     h_matrix_dressed(j,i,istate) += dressing_column_h(j,istate)
+     h_matrix_dressed(i,j,dressing_state) += dressing_column_h(j,dressing_state)
+     h_matrix_dressed(j,i,dressing_state) += dressing_column_h(j,dressing_state)
    enddo
-   h_matrix_dressed(i,i,istate) -= dressing_column_h(i,istate)
+   h_matrix_dressed(i,i,dressing_state) -= dressing_column_h(i,dressing_state)
  enddo
 END_PROVIDER
 
diff --git a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
index 9d1c50d4..3d95d6b0 100644
--- a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
+++ b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
@@ -513,7 +513,7 @@ END_PROVIDER
   double precision :: norm_left, stato
   integer, external :: pt2_find  
 
-  pt2_weight(1) = psi_coef_generators(1,pt2_stoch_istate)**2
+  pt2_weight(1)  = psi_coef_generators(1,pt2_stoch_istate)**2
   pt2_cweight(1) = psi_coef_generators(1,pt2_stoch_istate)**2
   
   do i=1,N_det_generators
diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f
index acda9fa6..d4fae71a 100644
--- a/plugins/Full_CI_ZMQ/selection.irp.f
+++ b/plugins/Full_CI_ZMQ/selection.irp.f
@@ -5,7 +5,6 @@ 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
 
diff --git a/plugins/MRCC_Utils/mrcc_general.irp.f b/plugins/MRPT_Utils/set_as_holes_and_particles.irp.f
similarity index 100%
rename from plugins/MRCC_Utils/mrcc_general.irp.f
rename to plugins/MRPT_Utils/set_as_holes_and_particles.irp.f
diff --git a/plugins/Psiref_Utils/psi_ref_utils.irp.f b/plugins/Psiref_Utils/psi_ref_utils.irp.f
index 95c993f0..c59bbd9f 100644
--- a/plugins/Psiref_Utils/psi_ref_utils.irp.f
+++ b/plugins/Psiref_Utils/psi_ref_utils.irp.f
@@ -325,3 +325,17 @@ BEGIN_PROVIDER [double precision, ref_hamiltonian_matrix, (n_det_ref,n_det_ref)]
  enddo
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ integer, idx_non_ref_from_sorted, (N_det) ]
+  implicit none
+  integer :: i,inpsisor
+  
+  idx_non_ref_from_sorted = 0
+
+  do i=1,N_det
+    inpsisor = psi_det_sorted_order(i)
+    if(inpsisor <= 0) stop "idx_non_ref_from_sorted"
+    idx_non_ref_from_sorted(inpsisor) = idx_non_ref_rev(i)
+  end do
+END_PROVIDER
+
diff --git a/plugins/dress_zmq/alpha_factory.irp.f b/plugins/dress_zmq/alpha_factory.irp.f
index f1a3a8e9..114d4dbb 100644
--- a/plugins/dress_zmq/alpha_factory.irp.f
+++ b/plugins/dress_zmq/alpha_factory.irp.f
@@ -20,19 +20,16 @@ subroutine alpha_callback(delta_ij_loc, i_generator, subset,iproc)
 end subroutine
 
 
- BEGIN_PROVIDER [ integer, psi_from_sorted, (N_det) ]
-&BEGIN_PROVIDER [ integer, idx_non_ref_from_sorted, (N_det) ]
+BEGIN_PROVIDER [ integer, psi_from_sorted, (N_det) ]
   implicit none
   integer :: i,inpsisor
   
-  idx_non_ref_from_sorted = 0
   psi_from_sorted = 0
 
   do i=1,N_det
     psi_from_sorted(psi_det_sorted_order(i)) = i
     inpsisor = psi_det_sorted_order(i)
     if(inpsisor <= 0) stop "idx_non_ref_from_sorted"
-    idx_non_ref_from_sorted(inpsisor) = idx_non_ref_rev(i)
   end do
 END_PROVIDER
 
@@ -41,7 +38,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
   use bitmasks
   implicit none
   BEGIN_DOC
-!            WARNING /!\ : It is assumed that the generators and selectors are psi_det_sorted
+! TODO
   END_DOC
   
   double precision,intent(inout) :: delta_ij_loc(N_states,N_det,2) 
@@ -379,6 +376,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
             call alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexes, indexes_end, abuf, siz, iproc)
             
             !call dress_with_alpha_buffer(delta_ij_loc, minilist, interesting(0), abuf, n)
+
           end if
         enddo
         if(s1 /= s2) monoBdo = .false.
@@ -489,10 +487,7 @@ subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexe
       !APPLY PART
       if(st4 > 1) then
         call apply_particles(mask, s1, i, s2, j, alpha, ok, N_int)
-        !if(.not. ok) stop "non existing alpha......"
-        !print *, "willcall", st4-1, size(labuf)
-        call dress_with_alpha_buffer(delta_ij_loc, labuf, det_minilist, st4-1, alpha, iproc)
-        !call dress_with_alpha_buffer(delta_ij_loc, abuf, siz, alpha, 1)
+        call dress_with_alpha_buffer(N_states, N_det, N_int, delta_ij_loc, labuf, det_minilist, st4-1, alpha, iproc)
       end if
     end do
     
diff --git a/plugins/dress_zmq/dress_general.irp.f b/plugins/dress_zmq/dress_general.irp.f
index 29dc086d..0bf7e715 100644
--- a/plugins/dress_zmq/dress_general.irp.f
+++ b/plugins/dress_zmq/dress_general.irp.f
@@ -1,6 +1,6 @@
 
 
-subroutine run(N_st,energy)
+subroutine run_dressing(N_st,energy)
   implicit none
   
   integer, intent(in) :: N_st 
@@ -13,8 +13,6 @@ subroutine run(N_st,energy)
   
   integer :: n_it_dress_max
   double precision :: thresh_dress
-  double precision, allocatable :: lambda(:)
-  allocate (lambda(N_states))
 
   thresh_dress = thresh_dressed_ci
   n_it_dress_max = n_it_max_dressed_ci
@@ -33,7 +31,6 @@ subroutine run(N_st,energy)
     E_new = 0.d0
     delta_E = 1.d0
     iteration = 0
-    lambda = 1.d0
     do while (delta_E > thresh_dress)
       iteration += 1
       print *,  '===============================================' 
@@ -44,7 +41,7 @@ subroutine run(N_st,energy)
       do i=1,N_st
         call write_double(6,ci_energy_dressed(i),"Energy")
       enddo
-      call diagonalize_ci_dressed(lambda)
+      call diagonalize_ci_dressed
       E_new = dress_e0_denominator(1) !sum(ci_energy_dressed(1:N_states))
 
       delta_E = (E_new - E_old)/dble(N_states)
@@ -63,72 +60,3 @@ subroutine run(N_st,energy)
   energy(1:N_st) = ci_energy_dressed(1:N_st)
 end
 
-
-subroutine print_cas_coefs
-  implicit none
-
-  integer :: i,j
-  print *,  'CAS'
-  print *,  '==='
-  do i=1,N_det_cas
-    print *,  (psi_cas_coef(i,j), j=1,N_states)
-    call debug_det(psi_cas(1,1,i),N_int)
-  enddo
-  call write_double(6,ci_energy(1),"Initial CI energy")
-
-end
-
-
-subroutine run_pt2(N_st,energy,pt2) 
-  implicit none 
-  integer :: i,j,k 
-  integer, intent(in)          :: N_st 
-  double precision, intent(in) :: energy(N_st) 
-  double precision :: pt2(N_st)
-  double precision :: norm_pert(N_st),H_pert_diag(N_st)
-  
-  pt2 = 0d0
-  
-  print*,'Last iteration only to compute the PT2' 
-  
-  N_det_generators = N_det_cas
-  N_det_selectors = N_det_non_ref
-
-  do i=1,N_det_generators
-    do k=1,N_int
-      psi_det_generators(k,1,i) = psi_ref(k,1,i)
-      psi_det_generators(k,2,i) = psi_ref(k,2,i)
-    enddo
-    do k=1,N_st
-      psi_coef_generators(i,k) = psi_ref_coef(i,k)
-    enddo
-  enddo
-  do i=1,N_det
-    do k=1,N_int
-      psi_selectors(k,1,i) = psi_det_sorted(k,1,i)
-      psi_selectors(k,2,i) = psi_det_sorted(k,2,i)
-    enddo
-    do k=1,N_st
-      psi_selectors_coef(i,k) = psi_coef_sorted(i,k)
-    enddo
-  enddo
-
-  SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
-  SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized
-
-  call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag,  N_st) 
-  
-!  call ezfio_set_full_ci_energy_pt2(energy+pt2)
-
-  print *,  'Final step' 
-  print *,  'N_det    = ', N_det 
-  print *,  'N_states = ', N_states 
-  print *,  'PT2      = ', pt2 
-  print *,  'E        = ', energy 
-  print *,  'E+PT2    = ', energy+pt2 
-  print *,  '-----' 
-
-!  call ezfio_set_dress_zmq_energy_pt2(energy(1)+pt2(1))
-
-end 
-
diff --git a/plugins/dress_zmq/dress_stoch_routines.irp.f b/plugins/dress_zmq/dress_stoch_routines.irp.f
index 3c2877ab..c2557758 100644
--- a/plugins/dress_zmq/dress_stoch_routines.irp.f
+++ b/plugins/dress_zmq/dress_stoch_routines.irp.f
@@ -10,10 +10,10 @@ subroutine ZMQ_dress(E, dress, delta, delta_s2, relative_error)
   implicit none
   
   character(len=64000)           :: task
-
+  
   integer(ZMQ_PTR)               :: zmq_to_qp_run_socket, zmq_socket_pull
   integer, external              :: omp_get_thread_num
-  double precision, intent(in)   :: relative_error, E
+  double precision, intent(in)   :: E(N_states), relative_error
   double precision, intent(out)  :: dress(N_states)
   double precision, intent(out)  :: delta(N_states, N_det)
   double precision, intent(out)  :: delta_s2(N_states, N_det)
@@ -23,106 +23,112 @@ subroutine ZMQ_dress(E, dress, delta, delta_s2, relative_error)
   
   double precision, external     :: omp_get_wtime
   double precision               :: time
-  double precision               :: w(N_states)
-  integer, external :: add_task_to_taskserver
-
-
-  provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral dress_weight psi_selectors
-
-  !!!!!!!!!!!!!!! demander a TOTO !!!!!!!
-  w(:) = 0.d0
-  w(dress_stoch_istate) = 1.d0
-  !call update_psi_average_norm_contrib(w)
+  integer, external              :: add_task_to_taskserver
+  double precision               :: state_average_weight_save(N_states)
   
   
-  
-  
-  print *, '========== ================= ================= ================='
-  print *, ' Samples        Energy         Stat. Error         Seconds      '
-  print *, '========== ================= ================= ================='
-  
-
-  call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull, 'dress')
-  
-  integer, external              :: zmq_put_psi
-  integer, external              :: zmq_put_N_det_generators
-  integer, external              :: zmq_put_N_det_selectors
-  integer, external              :: zmq_put_dvector
-  integer, external              :: zmq_set_running 
-  if (zmq_put_psi(zmq_to_qp_run_socket,1) == -1) then
-    stop 'Unable to put psi on ZMQ server'
-  endif
-  if (zmq_put_N_det_generators(zmq_to_qp_run_socket, 1) == -1) then
-    stop 'Unable to put N_det_generators on ZMQ server'
-  endif
-  if (zmq_put_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) then
-    stop 'Unable to put N_det_selectors on ZMQ server'
-  endif
-  if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',dress_e0_denominator,size(dress_e0_denominator)) == -1) then
-    stop 'Unable to put energy on ZMQ server'
-  endif
-
-  integer(ZMQ_PTR), external     :: new_zmq_to_qp_run_socket
-  integer                        :: ipos
-  ipos=1
-  do i=1,N_dress_jobs
-    if(dress_jobs(i) > fragment_first) then
-      write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') 0, dress_jobs(i)
-      ipos += 20
-      if (ipos > 63980) then
-        if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
-          stop 'Unable to add task to task server'
-        endif
-
-        ipos=1
-      endif
-    else
-      do j=1,fragment_count
-        write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') j, dress_jobs(i)
+  state_average_weight_save(:) = state_average_weight(:)
+  do dress_stoch_istate=1,N_states
+    SOFT_TOUCH dress_stoch_istate
+    state_average_weight(:) = 0.d0
+    state_average_weight(dress_stoch_istate) = 1.d0
+    TOUCH state_average_weight
+    
+    provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral dress_weight psi_selectors
+    
+    
+    print *, '========== ================= ================= ================='
+    print *, ' Samples        Energy         Stat. Error         Seconds      '
+    print *, '========== ================= ================= ================='
+    
+    
+    call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull, 'dress')
+    
+    integer, external              :: zmq_put_psi
+    integer, external              :: zmq_put_N_det_generators
+    integer, external              :: zmq_put_N_det_selectors
+    integer, external              :: zmq_put_dvector
+    integer, external              :: zmq_set_running
+    if (zmq_put_psi(zmq_to_qp_run_socket,1) == -1) then
+      stop 'Unable to put psi on ZMQ server'
+    endif
+    if (zmq_put_N_det_generators(zmq_to_qp_run_socket, 1) == -1) then
+      stop 'Unable to put N_det_generators on ZMQ server'
+    endif
+    if (zmq_put_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) then
+      stop 'Unable to put N_det_selectors on ZMQ server'
+    endif
+    if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',dress_e0_denominator,size(dress_e0_denominator)) == -1) then
+      stop 'Unable to put energy on ZMQ server'
+    endif
+    
+    integer(ZMQ_PTR), external     :: new_zmq_to_qp_run_socket
+    integer                        :: ipos
+    ipos=1
+    do i=1,N_dress_jobs
+      if(dress_jobs(i) > fragment_first) then
+        write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') 0, dress_jobs(i)
         ipos += 20
         if (ipos > 63980) then
           if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
             stop 'Unable to add task to task server'
           endif
+          
           ipos=1
         endif
-      end do
-    end if
-  end do
-  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 to task server'
-        endif
-  endif
-      if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
-        print *,  irp_here, ': Failed in zmq_set_running'
+      else
+        do j=1,fragment_count
+          write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') j, dress_jobs(i)
+          ipos += 20
+          if (ipos > 63980) then
+            if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
+              stop 'Unable to add task to task server'
+            endif
+            ipos=1
+          endif
+        end do
+      end if
+    end do
+    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 to task server'
       endif
-
-  !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc+1)                &
-      !$OMP  PRIVATE(i)
-  i = omp_get_thread_num()
-  if (i==0) then
-    call dress_collector(zmq_socket_pull,E, relative_error, delta, delta_s2, dress)
-  else
-    call dress_slave_inproc(i)
-  endif
-  !$OMP END PARALLEL
-  call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'dress')
+    endif
+    if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
+      print *,  irp_here, ': Failed in zmq_set_running'
+    endif
+    
+    !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc+1)              &
+        !$OMP  PRIVATE(i)
+    i = omp_get_thread_num()
+    if (i==0) then
+      call dress_collector(zmq_socket_pull,E, relative_error, delta, delta_s2, dress,&
+          dress_stoch_istate)
+    else
+      call dress_slave_inproc(i)
+    endif
+    !$OMP END PARALLEL
+    call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'dress')
+    
+    print *, '========== ================= ================= ================='
+  enddo
+  FREE dress_stoch_istate
+  state_average_weight(:) = state_average_weight_save(:)
+  TOUCH state_average_weight
   
-  print *, '========== ================= ================= ================='
 end subroutine
 
 
 subroutine dress_slave_inproc(i)
   implicit none
   integer, intent(in)            :: i
-
+  
   call run_dress_slave(1,i,dress_e0_denominator)
 end
 
 
 
-subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, dress)
+subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, dress, istate)
   use f77_zmq
   use bitmasks
   implicit none
@@ -130,8 +136,9 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
   
 
   integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull
+  integer, intent(in)            :: istate
 
-  double precision, intent(in) :: relative_error, E
+  double precision, intent(in)   :: relative_error, E(N_states)
   double precision, intent(out)  :: dress(N_states)
   double precision, allocatable  :: cp(:,:,:,:)
 
@@ -197,9 +204,9 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
     
     !!!!! A VERIFIER !!!!!
     do i_state=1,N_states
-    do i=1, N_det
-      dress_mwen(i_state) += delta_loc(i_state, i, 1) * psi_coef(i, i_state)
-    end do
+      do i=1, N_det
+        dress_mwen(i_state) += delta_loc(i_state, i, 1) * psi_coef(i, i_state)
+      end do
     end do
       
     dress_detail(:, ind) += dress_mwen(:)
@@ -210,12 +217,8 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
         integer :: toothMwen
         logical :: fracted
         fac = cps(ind, j) / cps_N(j) * dress_weight_inv(ind) * comb_step
-        do k=1,N_det
-        do i_state=1,N_states
-          cp(i_state,k,j,1) += delta_loc(i_state,k,1) * fac
-          cp(i_state,k,j,2) += delta_loc(i_state,k,2) * fac
-        end do
-        end do
+        cp(1:N_states,1:N_det,j,1) += delta_loc(1:N_states,1:N_det,1) * fac
+        cp(1:N_states,1:N_det,j,2) += delta_loc(1:N_states,1:N_det,2) * fac
       end if
     end do
     toothMwen = tooth_of_det(ind)
@@ -223,13 +226,13 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
     if(fracted) fracted = (ind == first_det_of_teeth(toothMwen))
     
     if(fracted) then
-      delta_det(:,:,toothMwen-1, 1) += delta_loc(:,:,1) * (1d0-fractage(toothMwen))
-      delta_det(:,:,toothMwen-1, 2) += delta_loc(:,:,2) * (1d0-fractage(toothMwen))
-      delta_det(:,:,toothMwen, 1) += delta_loc(:,:,1) * (fractage(toothMwen))
-      delta_det(:,:,toothMwen, 2) += delta_loc(:,:,2) * (fractage(toothMwen))
+      delta_det(1:N_states,1:N_det,toothMwen-1, 1) = delta_det(1:N_states,1:N_det,toothMwen-1, 1) + delta_loc(1:N_states,1:N_det,1) * (1d0-fractage(toothMwen))
+      delta_det(1:N_states,1:N_det,toothMwen-1, 2) = delta_det(1:N_states,1:N_det,toothMwen-1, 2) + delta_loc(1:N_states,1:N_det,2) * (1d0-fractage(toothMwen))
+      delta_det(1:N_states,1:N_det,toothMwen  , 1) = delta_det(1:N_states,1:N_det,toothMwen  , 1) + delta_loc(1:N_states,1:N_det,1) * (fractage(toothMwen))
+      delta_det(1:N_states,1:N_det,toothMwen  , 2) = delta_det(1:N_states,1:N_det,toothMwen  , 2) + delta_loc(1:N_states,1:N_det,2) * (fractage(toothMwen))
     else
-      delta_det(:,:,toothMwen, 1) += delta_loc(:,:,1)
-      delta_det(:,:,toothMwen, 2) += delta_loc(:,:,2)
+      delta_det(1:N_states,1:N_det,toothMwen  , 1) = delta_det(1:N_states,1:N_det,toothMwen  , 1) + delta_loc(1:N_states,1:N_det,1)
+      delta_det(1:N_states,1:N_det,toothMwen  , 2) = delta_det(1:N_states,1:N_det,toothMwen  , 2) + delta_loc(1:N_states,1:N_det,2)
     end if
 
     parts_to_get(ind) -= 1
@@ -265,25 +268,24 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
       double precision :: su, su2, eqt, avg, E0, val
       integer, external :: zmq_abort
 
-      su = 0d0
+      su  = 0d0
       su2 = 0d0
       
       do i=1, int(cps_N(cur_cp))
-        call get_comb_val(comb(i), dress_detail, cur_cp, val)
-        su += val
-        su2 += val**2
+        call get_comb_val(comb(i), dress_detail, cur_cp, val, istate)
+        su  += val
+        su2 += val*val
       end do
       avg = su / cps_N(cur_cp)
-      eqt = dsqrt( ((su2 / cps_N(cur_cp)) - avg**2) / cps_N(cur_cp) )
-      E0 = sum(dress_detail(1, :first_det_of_teeth(cp_first_tooth(cur_cp))-1))
+      eqt = dsqrt( ((su2 / cps_N(cur_cp)) - avg*avg) / cps_N(cur_cp) )
+      E0 = sum(dress_detail(istate, :first_det_of_teeth(cp_first_tooth(cur_cp))-1))
       if(cp_first_tooth(cur_cp) <= comb_teeth) then
-        E0 = E0 + dress_detail(1, first_det_of_teeth(cp_first_tooth(cur_cp))) * (1d0-fractage(cp_first_tooth(cur_cp)))
+        E0 = E0 + dress_detail(istate, first_det_of_teeth(cp_first_tooth(cur_cp))) * (1d0-fractage(cp_first_tooth(cur_cp)))
       end if
       call wall_time(time)
-      if ((dabs(eqt) < relative_error .and. cps_N(cur_cp) >= 30)  .or. total_computed == N_det_generators) then
+      if ((dabs(eqt) < relative_error .and. cps_N(cur_cp) >= 5)  .or. total_computed == N_det_generators) then
         ! Termination
-        !print '(G10.3, 2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, ''
-!        print *, "GREPME", cur_cp, E+E0+avg, eqt, time-time0, total_computed
+        print '(2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', avg+E(istate)+E0, eqt, time-time0, ''
         if (zmq_abort(zmq_to_qp_run_socket) == -1) then
           call sleep(1)
           if (zmq_abort(zmq_to_qp_run_socket) == -1) then
@@ -293,31 +295,30 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
       else
         if (cur_cp > old_cur_cp) then
           old_cur_cp = cur_cp
-!          print *, "GREPME", cur_cp, E+E0+avg, eqt, time-time0, total_computed
-          !print '(G10.3, 2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, ''
+          print '(2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', avg+E(istate)+E0, eqt, time-time0, ''
         endif
       endif
     end if
   end do pullLoop
 
   if(total_computed == N_det_generators) then
-    delta = 0d0
-    delta_s2 = 0d0
+    delta   (1:N_states,1:N_det) = 0d0
+    delta_s2(1:N_states,1:N_det) = 0d0
     do i=comb_teeth+1,0,-1
-      delta += delta_det(:,:,i,1)
-      delta_s2 += delta_det(:,:,i,2)
+      delta   (1:N_states,1:N_det) = delta   (1:N_states,1:N_det) + delta_det(1:N_states,1:N_det,i,1)
+      delta_s2(1:N_states,1:N_det) = delta_s2(1:N_states,1:N_det) + delta_det(1:N_states,1:N_det,i,2)
     end do
   else
 
-    delta = cp(:,:,cur_cp,1)
-    delta_s2 = cp(:,:,cur_cp,2)
+    delta   (1:N_states,1:N_det) = cp(1:N_states,1:N_det,cur_cp,1)
+    delta_s2(1:N_states,1:N_det) = cp(1:N_states,1:N_det,cur_cp,2)
     do i=cp_first_tooth(cur_cp)-1,0,-1
-      delta += delta_det(:,:,i,1)
-      delta_s2 += delta_det(:,:,i,2)
+      delta   (1:N_states,1:N_det) = delta   (1:N_states,1:N_det) + delta_det(1:N_states,1:N_det,i,1)
+      delta_s2(1:N_states,1:N_det) = delta_s2(1:N_states,1:N_det) + delta_det(1:N_states,1:N_det,i,2)
     end do
 
   end if
-  dress(1) = E
+  dress(istate) = E(istate)+E0
   call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
 end subroutine
 
@@ -353,10 +354,14 @@ end function
 &BEGIN_PROVIDER [ integer, comb_teeth ]
 &BEGIN_PROVIDER [ integer, N_cps_max ]
   implicit none
-  comb_teeth = 16
-  N_cps_max = 32
+  BEGIN_DOC
+! N_cps_max : max number of checkpoints
+!
+! gen_per_cp : number of generators per checkpoint
+  END_DOC
+  comb_teeth = 64
+  N_cps_max = 64
   gen_per_cp = (N_det_generators / N_cps_max) + 1
-  N_cps_max += 1
 END_PROVIDER
 
 
@@ -457,9 +462,9 @@ END_PROVIDER
 END_PROVIDER
 
 
-subroutine get_comb_val(stato, detail, cur_cp, val)
+subroutine get_comb_val(stato, detail, cur_cp, val, istate)
   implicit none
-  integer, intent(in)           :: cur_cp
+  integer, intent(in)           :: cur_cp, istate
   integer                       :: first
   double precision, intent(in)  :: stato, detail(N_states, N_det_generators)
   double precision, intent(out) :: val
@@ -475,9 +480,9 @@ subroutine get_comb_val(stato, detail, cur_cp, val)
     !DIR$ FORCEINLINE
     k = dress_find(curs, dress_cweight,size(dress_cweight), first_det_of_teeth(j), first_det_of_teeth(j+1))
     if(k == first_det_of_teeth(first)) then
-     val += detail(1, k) * dress_weight_inv(k) * comb_step * fractage(first)
+     val += detail(istate, k) * dress_weight_inv(k) * comb_step * fractage(first)
     else
-     val += detail(1, k) * dress_weight_inv(k) * comb_step
+     val += detail(istate, k) * dress_weight_inv(k) * comb_step
     end if
 
     curs -= comb_step
@@ -528,10 +533,10 @@ subroutine add_comb(com, computed, cp, N, tbc)
 end subroutine
 
 
- BEGIN_PROVIDER [ integer, dress_stoch_istate ]
-   implicit none
-   dress_stoch_istate = 1
- END_PROVIDER
+BEGIN_PROVIDER [ integer, dress_stoch_istate ]
+  implicit none
+  dress_stoch_istate = 1
+END_PROVIDER
 
  
  BEGIN_PROVIDER [ double precision, dress_weight, (N_det_generators) ]
@@ -548,7 +553,7 @@ end subroutine
   double precision :: norm_left, stato
   integer, external :: dress_find  
 
-  dress_weight(1) = psi_coef_generators(1,dress_stoch_istate)**2
+  dress_weight(1)  = psi_coef_generators(1,dress_stoch_istate)**2
   dress_cweight(1) = psi_coef_generators(1,dress_stoch_istate)**2
   
   do i=1,N_det_generators
@@ -575,8 +580,8 @@ end subroutine
   
   comb_step = 1d0/dfloat(comb_teeth)
   first_det_of_comb = 1
-  do i=1,N_det_generators
-    if(dress_weight(i)/norm_left < .25d0*comb_step) then
+  do i=1,min(100,N_det_generators)
+    if(dress_weight(i)/norm_left < comb_step) then
       first_det_of_comb = i
       exit
     end if
diff --git a/plugins/dress_zmq/dress_zmq_routines.irp.f b/plugins/dress_zmq/dress_zmq_routines.irp.f
index a18bc882..4dc75236 100644
--- a/plugins/dress_zmq/dress_zmq_routines.irp.f
+++ b/plugins/dress_zmq/dress_zmq_routines.irp.f
@@ -15,10 +15,7 @@ subroutine dress_zmq()
     enddo
     SOFT_TOUCH psi_coef 
   endif
-  call run(N_states,energy)
-  if(do_pt2)then
-    call run_pt2(N_states,energy)
-  endif
+  call run_dressing(N_states,energy)
   deallocate(energy)
 end
 
diff --git a/plugins/dress_zmq/dressing.irp.f b/plugins/dress_zmq/dressing.irp.f
index 1736cc76..f835ffbf 100644
--- a/plugins/dress_zmq/dressing.irp.f
+++ b/plugins/dress_zmq/dressing.irp.f
@@ -5,9 +5,9 @@ BEGIN_PROVIDER [ integer, N_dress_teeth ]
   N_dress_teeth = 10
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, dress_norm_acc, (0:N_det_non_ref, N_states) ]
-&BEGIN_PROVIDER [ double precision, dress_norm, (0:N_det_non_ref, N_states) ]
-&BEGIN_PROVIDER [ double precision, dress_teeth_size, (0:N_det_non_ref, N_states) ]
+BEGIN_PROVIDER [ double precision, dress_norm_acc, (0:N_det, N_states) ]
+&BEGIN_PROVIDER [ double precision, dress_norm, (0:N_det, N_states) ]
+&BEGIN_PROVIDER [ double precision, dress_teeth_size, (0:N_det, N_states) ]
 &BEGIN_PROVIDER [ integer, dress_teeth, (0:N_dress_teeth+1, N_states) ]
   implicit none
   integer :: i, j, st, nt
@@ -43,11 +43,11 @@ BEGIN_PROVIDER [ double precision, dress_norm_acc, (0:N_det_non_ref, N_states) ]
       end if
     end do
     if(nt > N_dress_teeth+1) then
-      print *, "foireouse dress_teeth", nt, dress_teeth(nt,st), N_det_non_ref
+      print *, "foireouse dress_teeth", nt, dress_teeth(nt,st), N_det
       stop
     end if
 
-    dress_teeth(N_dress_teeth+1,st) = N_det_non_ref+1
+    dress_teeth(N_dress_teeth+1,st) = N_det+1
     norm_loc = 0d0
     do i=N_dress_teeth, 0, -1
       dress_teeth_size(i,st) = dress_norm_acc(dress_teeth(i+1,st)-1,st) - dress_norm_acc(dress_teeth(i,st)-1, st)
@@ -64,39 +64,37 @@ END_PROVIDER
 
 
 
- BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det,2) ]
+BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det,2) ]
   use bitmasks
   implicit none
 
   integer                        :: i,j,k
   
   double precision, allocatable  :: dress(:), del(:,:), del_s2(:,:)
-  double precision               :: E_CI_before, relative_error
-  double precision, save :: errr = 0d0
+  double precision               :: E_CI_before(N_states), relative_error
+!  double precision, save         :: errr = 0d0
 
   allocate(dress(N_states), del(N_states, N_det), del_s2(N_states, N_det))
 
   delta_ij = 0d0
 
-  E_CI_before = dress_E0_denominator(1) + nuclear_repulsion
+  E_CI_before(:) = dress_E0_denominator(:) + nuclear_repulsion
   threshold_selectors = 1.d0
   threshold_generators = 1d0 
-  if(errr /= 0d0) then
-    errr = errr / 2d0 
-  else
-    errr = 1d-4
-  end if
-  relative_error = errr * 0d0
-  print *, "RELATIVE ERROR", relative_error
+!  if(errr /= 0d0) then
+!    errr = errr / 2d0 
+!  else
+!    errr = 1d-4
+!  end if
+  relative_error = 1.d-3
+  call write_double(6,relative_error,"Convergence of the stochastic algorithm")
+
   call ZMQ_dress(E_CI_before, dress, del, del_s2, abs(relative_error))
-  
   delta_ij(:,:,1) = del(:,:)
-  !delta_ij_s2(:,:,1) = del_s2(:,:)
   delta_ij(:,:,2) = del_s2(:,:)
-  !do i=N_det,1,-1
-  !  delta_ii(dress_stoch_istate,1) -= delta_ij(dress_stoch_istate, i, 1) / psi_ref_coef(1,dress_stoch_istate) * psi_coef(i, dress_stoch_istate)
-  !  delta_ii_s2(dress_stoch_istate,1) -= delta_ij_s2(dress_stoch_istate, i, 1) / psi_ref_coef(1,dress_stoch_istate) * psi_coef(i, dress_stoch_istate)
-  !end do
+
+  deallocate(dress, del, del_s2)
+  
 END_PROVIDER
 
 
diff --git a/plugins/dress_zmq/dressing_vector.irp.f b/plugins/dress_zmq/dressing_vector.irp.f
index d3e465bd..494e7c4b 100644
--- a/plugins/dress_zmq/dressing_vector.irp.f
+++ b/plugins/dress_zmq/dressing_vector.irp.f
@@ -5,25 +5,28 @@
  BEGIN_DOC
  ! Null dressing vectors
  END_DOC
- dressing_column_h(:,:) = 0.d0
- dressing_column_s(:,:) = 0.d0
 
- integer :: i,ii,k,j,jj, l
+ integer :: i,ii,k,j, l
  double precision :: f, tmp
  double precision, external :: u_dot_v
  
+ dressing_column_h(:,:) = 0.d0
+ dressing_column_s(:,:) = 0.d0
+
  do k=1,N_states
    l = dressed_column_idx(k)
    f = 1.d0/psi_coef(l,k)
-   do jj = 1, n_det
-     j = jj !idx_non_ref(jj)
-     dressing_column_h(j,k) = delta_ij   (k,jj,1) * f
-     dressing_column_s(j,k) = delta_ij   (k,jj,2) * f!delta_ij_s2(k,jj)
+   do j = 1, n_det
+     dressing_column_h(j,k) = delta_ij(k,j,1) * f
+     dressing_column_s(j,k) = delta_ij(k,j,2) * f
    enddo
-   tmp = u_dot_v(dressing_column_h(1,k), psi_coef(1,k), N_det)
+   tmp = u_dot_v(dressing_column_h(1,k), psi_coef(1,k), N_det) &
+     - dressing_column_h(l,k) * psi_coef(l,k)
    dressing_column_h(l,k) -= tmp * f
-   tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det)
+   tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) &
+     - dressing_column_s(l,k) * psi_coef(l,k)
    dressing_column_s(l,k) -= tmp * f
  enddo
+
 END_PROVIDER
 
diff --git a/plugins/dress_zmq/energy.irp.f b/plugins/dress_zmq/energy.irp.f
index 0ab170f1..b8948219 100644
--- a/plugins/dress_zmq/energy.irp.f
+++ b/plugins/dress_zmq/energy.irp.f
@@ -11,10 +11,12 @@ BEGIN_PROVIDER [ double precision, dress_E0_denominator, (N_states) ]
  BEGIN_DOC
  ! E0 in the denominator of the dress
  END_DOC
+ integer :: i
  if (initialize_dress_E0_denominator) then
-  dress_E0_denominator(1:N_states) = psi_energy(1:N_states)
-! dress_E0_denominator(1:N_states) = HF_energy - nuclear_repulsion
-! dress_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
+  call u_0_H_u_0(dress_E0_denominator,psi_coef,N_det,psi_det,N_int,N_states,size(psi_coef,1))
+  do i=N_det+1,N_states
+    dress_E0_denominator(i) = 0.d0
+  enddo
   call write_double(6,dress_E0_denominator(1)+nuclear_repulsion, 'dress Energy denominator')
  else
    dress_E0_denominator = -huge(1.d0)
diff --git a/plugins/dress_zmq/run_dress_slave.irp.f b/plugins/dress_zmq/run_dress_slave.irp.f
index 08d8af3d..b0896c00 100644
--- a/plugins/dress_zmq/run_dress_slave.irp.f
+++ b/plugins/dress_zmq/run_dress_slave.irp.f
@@ -46,7 +46,7 @@ subroutine run_dress_slave(thread,iproc,energy)
     return
   end if
   do i=1,N_states
-    div(i) = psi_ref_coef(dressed_column_idx(i), i)
+    div(i) = psi_coef(dressed_column_idx(i), i)
   end do
   do
     call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task)
@@ -56,14 +56,6 @@ subroutine run_dress_slave(thread,iproc,energy)
       delta_ij_loc = 0d0
       call alpha_callback(delta_ij_loc, i_generator, subset, iproc)
 
-      !!! SET DRESSING COLUMN?
-      !do i=1,N_det
-      !  do i_state=1,N_states
-      !    delta_ij_loc(i_state,i,1) = delta_ij_loc(i_state,i,1) / div(i_state) 
-      !    delta_ij_loc(i_state,i,2) = delta_ij_loc(i_state,i,2) / div(i_state)
-      !  end do
-      !end do
-
       call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
       call push_dress_results(zmq_socket_push, i_generator, delta_ij_loc, task_id)
     else
diff --git a/plugins/mrcc_sto/mrcc_sto.irp.f b/plugins/mrcc_sto/mrcc_sto.irp.f
index b8392f7f..cb92e3fa 100644
--- a/plugins/mrcc_sto/mrcc_sto.irp.f
+++ b/plugins/mrcc_sto/mrcc_sto.irp.f
@@ -16,7 +16,7 @@ end
 &BEGIN_PROVIDER [ integer, excs_ , (0:2,2,2,N_det,Nproc) ]
 &BEGIN_PROVIDER [ double precision, phases_, (N_det, Nproc) ]
 BEGIN_DOC
- ! temporay arrays for dress_with_alpha_buffer. Avoids realocation.
+ ! temporay arrays for dress_with_alpha_buffer. Avoids reallocation.
 END_DOC
 END_PROVIDER
 
diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f
index dc020fc9..3e5610c8 100644
--- a/src/Davidson/u0Hu0.irp.f
+++ b/src/Davidson/u0Hu0.irp.f
@@ -483,8 +483,11 @@ subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint,N_st,sze)
     call H_S2_u_0_nstates_zmq(v_0,s_0,u_1,N_states_diag,sze)
     deallocate(u_1)
   else
-    allocate (v_0(sze,N_st),s_0(sze,N_st))
-    call H_S2_u_0_nstates_openmp(v_0,s_0,u_0,N_st,sze)
+    allocate (v_0(n,N_st),s_0(n,N_st),u_1(n,N_st))
+    u_1(1:n,:) = u_0(1:n,:)
+    call H_S2_u_0_nstates_openmp(v_0,s_0,u_1,N_st,n)
+    u_0(1:n,:) = u_1(1:n,:)
+    deallocate(u_1)
   endif
   double precision :: norm
   do i=1,N_st