diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f
index de59e725..620605ff 100644
--- a/plugins/mrcepa0/dressing.irp.f
+++ b/plugins/mrcepa0/dressing.irp.f
@@ -1,6 +1,272 @@
 use bitmasks
 
 
+BEGIN_PROVIDER [ integer, N_mrcc_teeth ]
+  N_mrcc_teeth = 10
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, mrcc_norm_acc, (0:N_det_non_ref, N_states) ]
+&BEGIN_PROVIDER [ double precision, mrcc_norm, (0:N_det_non_ref, N_states) ]
+&BEGIN_PROVIDER [ double precision, mrcc_teeth_size, (0:N_det_non_ref, N_states) ]
+&BEGIN_PROVIDER [ integer, mrcc_teeth, (0:N_mrcc_teeth+1, N_states) ]
+  implicit none
+  integer :: i, j, st, nt
+  double precision :: norm_sto, jump, norm_mwen, norm_loc
+  
+  if(N_states /= 1) stop "mrcc_sto may not work with N_states /= 1"
+  
+  do st=1,N_states
+
+    !norm_non_ref = 1d0
+    !do i=1,N_det_ref
+    !  norm_non_ref -= psi_ref_coef(i,st)**2
+    !end do
+    
+    mrcc_teeth(0,st) = 1
+    !norm_non_sto = norm_non_ref
+    norm_sto = 1d0
+    do i=1,N_det
+      mrcc_teeth(1,st) = i
+      jump = (1d0 / dfloat(N_mrcc_teeth)) * norm_sto
+      if(psi_coef_generators(i,1)**2 < jump / 2d0) exit
+      !if(i==80) exit
+      norm_sto -= psi_coef_generators(i,1)**2
+    end do
+    
+    print *, "FIRST", mrcc_teeth(1,st)
+
+    norm_loc = 0d0
+    mrcc_norm_acc(0,st) = 0d0
+    nt = 1
+    
+    do i=1,mrcc_teeth(1,st)-1
+      mrcc_norm_acc(i,st) = mrcc_norm_acc(i-1,st) + psi_coef_generators(i,st)**2
+    end do
+
+    do i=mrcc_teeth(1,st), N_det_generators!-mrcc_teeth(1,st)+1
+      norm_mwen = psi_coef_generators(i,st)**2!-1+mrcc_teeth(1,st),st)**2
+      mrcc_norm_acc(i,st) = mrcc_norm_acc(i-1,st) + norm_mwen 
+      norm_loc += norm_mwen
+      if(norm_loc > (jump*dfloat(nt))) then
+        nt = nt + 1
+        mrcc_teeth(nt,st) = i
+      end if
+    end do
+    if(nt > N_mrcc_teeth+1) then
+      print *, "foireouse mrcc_teeth", nt, mrcc_teeth(nt,st), N_det_non_ref
+      stop
+    end if
+
+    mrcc_teeth(N_mrcc_teeth+1,st) = N_det_non_ref+1
+    !mrcc_norm_acc(:,st) = mrcc_norm_acc(:,st) / norm_non_ref
+    norm_loc = 0d0
+    do i=N_mrcc_teeth, 0, -1
+      mrcc_teeth_size(i,st) = mrcc_norm_acc(mrcc_teeth(i+1,st)-1,st) - mrcc_norm_acc(mrcc_teeth(i,st)-1, st)
+      mrcc_norm_acc(mrcc_teeth(i,st):mrcc_teeth(i+1,st)-1,st) -= mrcc_norm_acc(mrcc_teeth(i,st)-1, st)
+      mrcc_norm_acc(mrcc_teeth(i,st):mrcc_teeth(i+1,st)-1,st) = &
+          mrcc_norm_acc(mrcc_teeth(i,st):mrcc_teeth(i+1,st)-1,st) / mrcc_teeth_size(i,st)
+      mrcc_norm(mrcc_teeth(i,st), st) = mrcc_norm_acc(mrcc_teeth(i,st), st)
+      do j=mrcc_teeth(i,st)+1, mrcc_teeth(i+1,1)-1
+        mrcc_norm(j,1) = mrcc_norm_acc(j, st) - mrcc_norm_acc(j-1, st)
+      end do
+    end do
+  end do
+END_PROVIDER
+ 
+
+ BEGIN_PROVIDER [ double precision, delta_ij_mrcc_sto,(N_states,N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ii_mrcc_sto, (N_states, N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_sto, (N_states,N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ii_s2_mrcc_sto, (N_states, N_det_ref) ]
+  use bitmasks
+  implicit none
+  integer :: gen, h, p, n, t, i, j, h1, h2, p1, p2, s1, s2, iproc
+  integer(bit_kind) :: mask(N_int, 2), omask(N_int, 2)
+  integer(bit_kind),allocatable :: buf(:,:,:)
+  logical :: ok
+  logical, external :: detEq
+  integer, external :: omp_get_thread_num
+  double precision :: coefs(N_det_non_ref), myCoef
+  integer :: n_in_teeth
+  double precision :: contrib(N_states), curn, in_teeth_step, curlim, curnorm
+  
+  contrib = 0d0
+  read(*,*) n_in_teeth
+  !n_in_teeth = 2
+  in_teeth_step = 1d0 / dfloat(n_in_teeth)
+!double precision :: delta_ij_mrcc_tmp,(N_states,N_det_non_ref,N_det_ref) ]
+  !double precision :: delta_ii_mrcc_tmp, (N_states,N_det_ref) ]
+  !double precision :: delta_ij_s2_mrcc_tmp(N_states,N_det_non_ref,N_det_ref)
+  !double precision :: delta_ii_s2_mrcc_tmp(N_states, N_det_ref)
+
+  coefs = 0d0
+  coefs(:mrcc_teeth(1,1)-1) = 1d0
+  
+  do i=1,N_mrcc_teeth
+    print *, "TEETH SIZE", i, mrcc_teeth(i+1,1)-mrcc_teeth(i,1)
+    if(mrcc_teeth(i+1,1) - mrcc_teeth(i,1) <= n_in_teeth) then
+      coefs(mrcc_teeth(i,1):mrcc_teeth(i+1,1)-1) = 1d0
+    else if(.false.) then
+      curnorm = 0d0
+      curn = 0.5d0
+      curlim = curn / dfloat(n_in_teeth)
+      do j=mrcc_teeth(i,1), mrcc_teeth(i+1,1)-1
+        if(mrcc_norm_acc(j,1) >= curlim) then
+          coefs(j) = 1d0
+          curnorm += mrcc_norm(j,1)
+          do while(mrcc_norm_acc(j,1) > curlim)
+            curn += 1d0
+            curlim = curn / dfloat(n_in_teeth)
+          end do
+        end if
+      end do
+      do j=mrcc_teeth(i,1), mrcc_teeth(i+1,1)-1
+        coefs(j) = coefs(j) / curnorm ! 1d0 / norm computed in teeth
+      end do
+    else if(.true.) then
+      coefs(mrcc_teeth(i,1):mrcc_teeth(i,1)+n_in_teeth-1) = 1d0 / mrcc_norm_acc(mrcc_teeth(i,1)+n_in_teeth-1, 1)
+    else
+      curnorm = 0d0
+      n = mrcc_teeth(i+1,1) - mrcc_teeth(i,1)
+      do j=1,n_in_teeth
+        t = int((dfloat(j)-0.5d0) * dfloat(n) / dfloat(n_in_teeth)) + 1 + mrcc_teeth(i,1) - 1
+        curnorm += mrcc_norm(t,1)
+        coefs(t) = 1d0
+      end do
+      do j=mrcc_teeth(i,1), mrcc_teeth(i+1,1)-1
+        coefs(j) = coefs(j) / curnorm ! 1d0 / norm computed in teeth
+      end do
+    end if
+    !coefs(mrcc_teeth(i,1)) = 
+  end do
+
+  !coefs = coefs * dfloat(N_det_generators)
+
+
+  delta_ij_mrcc_sto = 0d0
+  delta_ii_mrcc_sto = 0d0
+  delta_ij_s2_mrcc_sto = 0d0
+  delta_ii_s2_mrcc_sto = 0d0
+  PROVIDE dij
+  provide hh_shortcut psi_det_size! lambda_mrcc
+  !$OMP PARALLEL DO default(none)  schedule(dynamic) &
+  !$OMP shared(psi_ref, psi_non_ref, hh_exists, pp_exists, N_int, hh_shortcut) &
+  !$OMP shared(N_det_generators, coefs,N_det_non_ref, N_det_ref, delta_ii_mrcc_sto, delta_ij_mrcc_sto) &
+  !$OMP shared(contrib,psi_det_generators, delta_ii_s2_mrcc_sto, delta_ij_s2_mrcc_sto) &
+  !$OMP private(i,j,curnorm,myCoef, h, n, mask, omask, buf, ok, iproc) 
+  do gen= 1,N_det_generators
+    if(coefs(gen) == 0d0) cycle
+    myCoef = coefs(gen)
+    allocate(buf(N_int, 2, N_det_non_ref))
+    iproc = omp_get_thread_num() + 1
+    if(mod(gen, 1000) == 0) print *, "mrcc_sto ", gen, "/", N_det_generators
+    
+    do h=1, hh_shortcut(0)
+      call apply_hole_local(psi_det_generators(1,1,gen), hh_exists(1, h), mask, ok, N_int)
+      if(.not. ok) cycle
+      omask = 0_bit_kind
+      if(hh_exists(1, h) /= 0) omask = mask
+      n = 1
+      do p=hh_shortcut(h), hh_shortcut(h+1)-1
+        call apply_particle_local(mask, pp_exists(1, p), buf(1,1,n), ok, N_int)
+        if(ok) n = n + 1
+        if(n > N_det_non_ref) stop "Buffer too small in MRCC..."
+      end do
+      n = n - 1
+      if(n /= 0) then
+        call mrcc_part_dress(delta_ij_mrcc_sto, delta_ii_mrcc_sto, delta_ij_s2_mrcc_sto, &
+            delta_ii_s2_mrcc_sto, gen,n,buf,N_int,omask,myCoef,contrib)
+      endif
+    end do
+    deallocate(buf)
+    end do
+    !$OMP END PARALLEL DO
+
+
+
+    curnorm = 0d0
+    do i=1,N_det_ref
+    do j=1,N_det_non_ref
+      curnorm += delta_ij_mrcc_sto(1, j, i)**2
+    end do
+    end do
+    print *, "NORM DELTA ", curnorm**0.5d0
+
+END_PROVIDER
+
+
+
+ BEGIN_PROVIDER [ double precision, delta_ij_cancel, (N_states,N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ii_cancel, (N_states, N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ij_s2_cancel, (N_states,N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ii_s2_cancel, (N_states, N_det_ref) ]
+  use bitmasks
+  implicit none
+
+
+  integer                         :: i_state, i, i_I, J, k, k2, k1, kk, ll, m,l, deg, ni, m2
+  integer                         :: n(2)
+  integer                         :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s, kn
+  logical                         :: ok
+  double precision                :: phase_ia, phase_Ik, phase_Jl, phase_Ji,phase_la, phase_ka, phase_tmp
+  double precision                :: Hka, Hla, Ska, Sla, tmp
+  double precision                :: diI, hIi, hJi, delta_JI, dkI, HkI,ci_inv(N_states), cj_inv(N_states)
+  double precision                :: contrib, contrib_s2, wall, iwall
+  integer, dimension(0:2,2,2)     :: exc_iI, exc_Ik, exc_IJ, exc
+  integer(bit_kind)               :: det_tmp(N_int, 2), det_tmp2(N_int, 2),inac, virt
+  integer, external               :: get_index_in_psi_det_sorted_bit, searchDet,detCmp
+  logical, external               :: is_in_wavefunction
+     
+  provide dij
+  
+  delta_ij_cancel = 0d0
+  delta_ii_cancel = 0d0
+
+  do i=1,N_det_ref
+    !$OMP PARALLEL DO default(shared) private(kk, k, blok, exc_Ik,det_tmp2,ok,deg,phase_Ik, l,ll) &
+    !$OMP private(contrib, contrib_s2, i_state)
+    do kk = 1, nlink(i)
+      k = det_cepa0_idx(linked(kk, i))
+      blok = blokMwen(kk, i)
+      call get_excitation(psi_ref(1,1,i),psi_non_ref(1,1,k),exc_Ik,deg,phase_Ik,N_int)
+      
+      do j=1,N_det_ref
+        if(j == i) cycle
+        call apply_excitation(psi_ref(1,1,J),exc_Ik,det_tmp2,ok,N_int)
+        if(.not. ok) cycle
+          
+        l = searchDet(det_cepa0(1,1,cepa0_shortcut(blok)), det_tmp2,cepa0_shortcut(blok+1)-cepa0_shortcut(blok), N_int)
+        if(l == -1) cycle
+        ll = cepa0_shortcut(blok)-1+l
+        l = det_cepa0_idx(ll)
+        ll = child_num(ll, J)
+     
+        do i_state = 1, N_states 
+          contrib =  (dij(j, l, i_state) - dij(i, k, i_state)) * delta_cas(i,j,i_state)! * Hla *phase_ia * phase_ik
+          contrib_s2 =  dij(j, l, i_state) - dij(i, k, i_state)! * Sla*phase_ia * phase_ik
+          if(dabs(psi_ref_coef(i,i_state)).ge.1.d-3) then
+            !$OMP ATOMIC
+            delta_ij_cancel(i_state,l,i) += contrib
+            !$OMP ATOMIC
+            delta_ij_s2_cancel(i_state,l,i) += contrib_s2
+            !$OMP ATOMIC
+            delta_ii_cancel(i_state,i) -= contrib / psi_ref_coef(i, i_state) * psi_non_ref_coef(l,i_state)
+            !$OMP ATOMIC
+            delta_ii_s2_cancel(i_state,i) -= contrib_s2 / psi_ref_coef(i, i_state) * psi_non_ref_coef(l,i_state)
+          else
+           !$OMP ATOMIC
+           delta_ij_cancel(i_state,l,i) += contrib * 0.5d0
+           !$OMP ATOMIC
+           delta_ij_s2_cancel(i_state,l,i) += contrib_s2 * 0.5d0
+         endif
+        end do
+      end do
+    end do
+    !$OMP END PARALLEL DO
+  enddo
+
+END_PROVIDER
+
 
  BEGIN_PROVIDER [ double precision, delta_ij_mrcc, (N_states,N_det_non_ref,N_det_ref) ]
 &BEGIN_PROVIDER [ double precision, delta_ii_mrcc, (N_states, N_det_ref) ]
@@ -14,15 +280,19 @@ use bitmasks
   logical :: ok
   logical, external :: detEq
   integer, external :: omp_get_thread_num
+  double precision :: contrib(N_states)
   
+  provide dij hh_shortcut psi_det_size  
+  
+  contrib = 0d0
   delta_ij_mrcc = 0d0
   delta_ii_mrcc = 0d0
   delta_ij_s2_mrcc = 0d0
   delta_ii_s2_mrcc = 0d0
-  PROVIDE dij
-  provide hh_shortcut psi_det_size! lambda_mrcc
+
+
   !$OMP PARALLEL DO default(none)  schedule(dynamic) &
-  !$OMP shared(psi_det_generators, N_det_generators, hh_exists, pp_exists, N_int, hh_shortcut) &
+  !$OMP shared(contrib,psi_det_generators, N_det_generators, hh_exists, pp_exists, N_int, hh_shortcut) &
   !$OMP shared(N_det_non_ref, N_det_ref, delta_ii_mrcc, delta_ij_mrcc, delta_ii_s2_mrcc, delta_ij_s2_mrcc) &
   !$OMP private(h, n, mask, omask, buf, ok, iproc) 
   do gen= 1, N_det_generators
@@ -43,7 +313,7 @@ use bitmasks
       n = n - 1
 
       if(n /= 0) then
-        call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc, delta_ij_s2_mrcc, delta_ii_s2_mrcc, gen,n,buf,N_int,omask)
+        call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc, delta_ij_s2_mrcc, delta_ii_s2_mrcc, gen,n,buf,N_int,omask,1d0,contrib)
       endif
 
     end do
@@ -59,7 +329,7 @@ END_PROVIDER
 ! end subroutine
 
 
-subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_generator,n_selected,det_buffer,Nint,key_mask)
+subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_generator,n_selected,det_buffer,Nint,key_mask,coef,contrib)
  use bitmasks
  implicit none
 
@@ -99,7 +369,10 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
   logical, external :: detEq, is_generable
   !double precision, external :: get_dij, get_dij_index
   double precision :: Delta_E_inv(N_states)
-  
+  double precision, intent(in) :: coef
+  double precision, intent(inout) :: contrib(N_states)
+  double precision :: sdress, hdress
+
   if (perturbative_triples) then
     PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
   endif
@@ -194,6 +467,288 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
       k_sd = idx_alpha(l_sd)
       call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd))
       call get_s2(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,sij_cache(k_sd))
+      !if(sij_cache(k_sd) /= 0D0) PRINT *, "SIJ ", sij_cache(k_sd) 
+    enddo
+
+    ! |I>
+    do i_I=1,N_det_ref
+      ! Find triples and quadruple grand parents
+      call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree1,Nint)
+      if (degree1 > 4) then
+        cycle
+      endif
+      
+      do i_state=1,N_states
+        dIa(i_state) = 0.d0
+      enddo
+      
+      ! <I|  <>  |alpha>
+      do k_sd=1,idx_alpha(0)
+
+        call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint)
+        if (degree > 2) then
+          cycle
+        endif
+
+        ! <I| /k\ |alpha>
+        
+        ! |l> = Exc(k -> alpha) |I>
+        call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree2,phase,Nint)
+        call decode_exc(exc,degree2,h1,p1,h2,p2,s1,s2)
+        do k=1,N_int
+          tmp_det(k,1) = psi_ref(k,1,i_I)
+          tmp_det(k,2) = psi_ref(k,2,i_I)
+        enddo
+        logical :: ok
+        call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint)
+        
+        do i_state=1,N_states
+          dIK(i_state) = dij(i_I, idx_alpha(k_sd), i_state)
+        enddo
+        
+        ! <I| \l/ |alpha>
+        do i_state=1,N_states
+          dka(i_state) = 0.d0
+        enddo
+
+        if (ok) then
+          do l_sd=k_sd+1,idx_alpha(0)
+            call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
+            if (degree == 0) then
+              call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
+              do i_state=1,N_states
+                dka(i_state) = dij(i_I, idx_alpha(l_sd), i_state) * phase * phase2
+              enddo
+              exit
+            endif
+          enddo
+
+        else if (perturbative_triples) then
+           ! Linked
+
+            hka = hij_cache(idx_alpha(k_sd))
+            if (dabs(hka) > 1.d-12) then
+              call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
+
+              do i_state=1,N_states
+                ASSERT (Delta_E_inv(i_state) < 0.d0)
+                dka(i_state) = hka / Delta_E_inv(i_state)
+              enddo
+            endif
+
+        endif
+
+        if (perturbative_triples.and. (degree2 == 1) ) then
+            call i_h_j(psi_ref(1,1,i_I),tmp_det,Nint,hka)
+            hka = hij_cache(idx_alpha(k_sd)) - hka
+            if (dabs(hka) > 1.d-12) then
+              call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv)
+              do i_state=1,N_states
+                ASSERT (Delta_E_inv(i_state) < 0.d0)
+                dka(i_state) = hka / Delta_E_inv(i_state)
+              enddo
+            endif
+
+        endif
+
+        do i_state=1,N_states
+          dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
+        enddo
+      enddo
+     
+      do i_state=1,N_states
+        ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state)
+      enddo
+      do l_sd=1,idx_alpha(0)
+        k_sd = idx_alpha(l_sd)
+        hla = hij_cache(k_sd)
+        sla = sij_cache(k_sd)
+        do i_state=1,N_states
+          dIa_hla(i_state,k_sd) = dIa(i_state) * hla * coef
+          dIa_sla(i_state,k_sd) = dIa(i_state) * sla * coef
+        enddo
+      enddo
+      do i_state=1,N_states
+        if(dabs(psi_ref_coef(1,i_state)).ge.1.d-3)then
+          do l_sd=1,idx_alpha(0)
+            k_sd = idx_alpha(l_sd)
+            p1 = 1
+            hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
+            sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
+            !$OMP ATOMIC
+            contrib(i_state) += hdress * psi_ref_coef(p1, i_state) * psi_non_ref_coef(k_sd, i_state)
+            !$OMP ATOMIC
+            delta_ij_(i_state,k_sd,p1) += hdress
+            !$OMP ATOMIC
+            !delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
+            delta_ii_(i_state,p1) -= hdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd)
+            !$OMP ATOMIC
+            delta_ij_s2_(i_state,k_sd,p1) += sdress
+            !$OMP ATOMIC
+            !delta_ii_s2_(i_state,i_I) = delta_ii_s2_(i_state,i_I) - dIa_sla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
+            delta_ii_s2_(i_state,p1) -= sdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd)
+          enddo
+        else
+          !stop "dress with coef < 1d-3"
+          delta_ii_(i_state,1)  = 0.d0
+          do l_sd=1,idx_alpha(0)
+            k_sd = idx_alpha(l_sd)
+            p1 = 1
+            hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
+            sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
+            !$OMP ATOMIC
+            delta_ij_(i_state,k_sd,p1) = delta_ij_(i_state,k_sd,p1) + 0.5d0*hdress
+            !$OMP ATOMIC
+            delta_ij_s2_(i_state,k_sd,p1) = delta_ij_s2_(i_state,k_sd,p1) + 0.5d0*sdress
+          enddo
+        endif
+      enddo
+    enddo
+  enddo
+  deallocate (dIa_hla,dIa_sla,hij_cache,sij_cache)
+  deallocate(miniList, idx_miniList)
+end
+
+
+
+subroutine mrcc_part_dress_1c(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_generator,n_selected,det_buffer,Nint,key_mask,contrib)
+ use bitmasks
+ implicit none
+
+  integer, intent(in)            :: i_generator,n_selected, Nint
+  double precision, intent(inout) :: delta_ij_(N_states,N_det_non_ref)
+  double precision, intent(inout) :: delta_ii_(N_states)
+  double precision, intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref)
+  double precision, intent(inout) :: delta_ii_s2_(N_states)
+
+  integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
+  integer                        :: i,j,k,l,m
+  integer,allocatable            :: idx_alpha(:), degree_alpha(:)
+  logical                        :: good, fullMatch
+
+  integer(bit_kind),allocatable  :: tq(:,:,:)
+  integer                        :: N_tq, c_ref ,degree1, degree2, degree
+
+  double precision               :: hIk, hla, hIl, sla, dIk(N_states), dka(N_states), dIa(N_states), hka
+  double precision, allocatable  :: dIa_hla(:,:), dIa_sla(:,:)
+  double precision               :: haj, phase, phase2
+  double precision               :: f(N_states), ci_inv(N_states)
+  integer                        :: exc(0:2,2,2)
+  integer                        :: h1,h2,p1,p2,s1,s2
+  integer(bit_kind)              :: tmp_det(Nint,2)
+  integer                        :: iint, ipos
+  integer                        :: i_state, k_sd, l_sd, i_I, i_alpha
+  
+  integer(bit_kind),allocatable  :: miniList(:,:,:)
+  integer(bit_kind),intent(in)   :: key_mask(Nint, 2)
+  integer,allocatable            :: idx_miniList(:)
+  integer                        :: N_miniList, ni, leng
+  double precision, allocatable  :: hij_cache(:), sij_cache(:)
+  
+  integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:)
+  integer, allocatable           :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:)
+  integer :: mobiles(2), smallerlist
+  logical, external :: detEq, is_generable
+  !double precision, external :: get_dij, get_dij_index
+  double precision :: Delta_E_inv(N_states)
+  double precision, intent(inout) :: contrib(N_states)
+  double precision :: sdress, hdress
+
+  if (perturbative_triples) then
+    PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
+  endif
+
+
+  leng = max(N_det_generators, N_det_non_ref)
+  allocate(miniList(Nint, 2, leng), tq(Nint,2,n_selected), idx_minilist(leng), hij_cache(N_det_non_ref), sij_cache(N_det_non_ref))
+  allocate(idx_alpha(0:psi_det_size), degree_alpha(psi_det_size))
+  call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
+  
+  allocate(ptr_microlist(0:mo_tot_num*2+1),  &
+      N_microlist(0:mo_tot_num*2) )
+  allocate(   microlist(Nint,2,N_minilist*4),               &
+    idx_microlist(N_minilist*4))
+  
+  if(key_mask(1,1) /= 0) then
+    call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
+    call filter_tq_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
+  else
+    call filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
+  end if
+  
+  
+  
+  deallocate(microlist, idx_microlist)
+  
+  allocate (dIa_hla(N_states,N_det_non_ref), dIa_sla(N_states,N_det_non_ref))
+  
+  ! |I>
+  
+  ! |alpha>
+  
+  if(N_tq > 0) then
+    call create_minilist(key_mask, psi_non_ref, miniList, idx_minilist, N_det_non_ref, N_minilist, Nint)
+    if(N_minilist == 0) return
+    
+    
+    if(sum(abs(key_mask(1:N_int,1))) /= 0) then
+      allocate(microlist_zero(Nint,2,N_minilist), idx_microlist_zero(N_minilist))
+      
+      allocate(   microlist(Nint,2,N_minilist*4),               &
+        idx_microlist(N_minilist*4))
+      call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
+      
+      
+      do i=0,mo_tot_num*2
+        do k=ptr_microlist(i),ptr_microlist(i+1)-1
+          idx_microlist(k) = idx_minilist(idx_microlist(k))
+        end do
+      end do
+      
+      do l=1,N_microlist(0)
+        do k=1,Nint
+          microlist_zero(k,1,l) = microlist(k,1,l)
+          microlist_zero(k,2,l) = microlist(k,2,l)
+        enddo
+        idx_microlist_zero(l) = idx_microlist(l)
+      enddo
+    end if
+  end if
+      
+ 
+  do i_alpha=1,N_tq
+   if(key_mask(1,1) /= 0) then
+      call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) 
+      
+      if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then
+        smallerlist = mobiles(1)
+      else
+        smallerlist = mobiles(2)
+      end if
+      
+    
+      do l=0,N_microlist(smallerlist)-1
+        microlist_zero(:,:,ptr_microlist(1) + l) = microlist(:,:,ptr_microlist(smallerlist) + l)
+        idx_microlist_zero(ptr_microlist(1) + l) = idx_microlist(ptr_microlist(smallerlist) + l)
+      end do
+      
+      call get_excitation_degree_vector(microlist_zero,tq(1,1,i_alpha),degree_alpha,Nint,N_microlist(smallerlist)+N_microlist(0),idx_alpha)
+      do j=1,idx_alpha(0)
+        idx_alpha(j) = idx_microlist_zero(idx_alpha(j))
+      end do
+      
+    else
+      call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
+      do j=1,idx_alpha(0)
+        idx_alpha(j) = idx_miniList(idx_alpha(j))
+      end do
+    end if
+    
+    do l_sd=1,idx_alpha(0)
+      k_sd = idx_alpha(l_sd)
+      call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd))
+      call get_s2(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,sij_cache(k_sd))
+      !if(sij_cache(k_sd) /= 0D0) PRINT *, "SIJ ", sij_cache(k_sd) 
     enddo
 
     ! |I>
@@ -295,26 +850,37 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
         enddo
       enddo
       do i_state=1,N_states
-        if(dabs(psi_ref_coef(i_I,i_state)).ge.1.d-3)then
+        if(dabs(psi_ref_coef(1,i_state)).ge.1.d-3)then
           do l_sd=1,idx_alpha(0)
             k_sd = idx_alpha(l_sd)
+            p1 = 1
+            hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
+            sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
             !$OMP ATOMIC
-            delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
+            contrib(i_state) += hdress * psi_ref_coef(p1, i_state) * psi_non_ref_coef(k_sd, i_state)
             !$OMP ATOMIC
-            delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
+            delta_ij_(i_state,k_sd) += hdress
             !$OMP ATOMIC
-            delta_ij_s2_(i_state,k_sd,i_I) = delta_ij_s2_(i_state,k_sd,i_I) + dIa_sla(i_state,k_sd)
+            !delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
+            delta_ii_(i_state) -= hdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd)
             !$OMP ATOMIC
-            delta_ii_s2_(i_state,i_I) = delta_ii_s2_(i_state,i_I) - dIa_sla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
+            delta_ij_s2_(i_state,k_sd) += sdress
+            !$OMP ATOMIC
+            !delta_ii_s2_(i_state,i_I) = delta_ii_s2_(i_state,i_I) - dIa_sla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
+            delta_ii_s2_(i_state) -= sdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd)
           enddo
         else
-          delta_ii_(i_state,i_I)  = 0.d0
+          !stop "dress with coef < 1d-3"
+          delta_ii_(i_state)  = 0.d0
           do l_sd=1,idx_alpha(0)
             k_sd = idx_alpha(l_sd)
+            p1 = 1
+            hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
+            sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state)
             !$OMP ATOMIC
-            delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + 0.5d0*dIa_hla(i_state,k_sd)
+            delta_ij_(i_state,k_sd) = delta_ij_(i_state,k_sd) + 0.5d0*hdress
             !$OMP ATOMIC
-            delta_ij_s2_(i_state,k_sd,i_I) = delta_ij_s2_(i_state,k_sd,i_I) + 0.5d0*dIa_sla(i_state,k_sd)
+            delta_ij_s2_(i_state,k_sd) = delta_ij_s2_(i_state,k_sd) + 0.5d0*sdress
           enddo
         endif
       enddo
@@ -325,6 +891,56 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
 end
 
 
+ BEGIN_PROVIDER [ double precision, mrcc_previous_E, (N_states) ]
+  implicit none
+  BEGIN_DOC
+    !energy difference between last two mrcc iterations
+  END_DOC
+  mrcc_previous_E(:) = ref_bitmask_energy
+ END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, delta_ij_mrcc_zmq, (N_states,N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ii_mrcc_zmq, (N_states, N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_zmq, (N_states,N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ii_s2_mrcc_zmq, (N_states, N_det_ref) ]  
+  use bitmasks
+  implicit none
+
+  integer                        :: i,j,k
+  
+  double precision, allocatable  :: mrcc(:)
+  double precision               :: E_CI_before, relative_error
+  double precision, save :: errr = 0d0
+
+  allocate(mrcc(N_states))
+
+
+  delta_ij_mrcc_zmq = 0d0
+  delta_ii_mrcc_zmq = 0d0
+  delta_ij_s2_mrcc_zmq = 0d0
+  delta_ii_s2_mrcc_zmq = 0d0
+
+  !call random_seed()
+  E_CI_before = mrcc_E0_denominator(1) + nuclear_repulsion
+  threshold_selectors = 1.d0
+  threshold_generators = 1d0 
+  !errr = errr / 2d0
+  if(errr /= 0d0) then
+    errr = errr / 2d0 ! (-mrcc_E0_denominator(1) + mrcc_previous_E(1)) / 1d1
+  else
+    errr = 1d-4
+  end if
+  relative_error = errr
+  print *, "RELATIVE ERROR", relative_error
+  call ZMQ_mrcc(E_CI_before, mrcc, delta_ij_mrcc_zmq, delta_ij_s2_mrcc_zmq, abs(relative_error))
+  !errr = 
+  mrcc_previous_E(:) = mrcc_E0_denominator(:)
+  do i=N_det_non_ref,1,-1
+    delta_ii_mrcc_zmq(:,1) -= delta_ij_mrcc_zmq(:, i, 1) / psi_ref_coef(1,1) * psi_non_ref_coef(i, 1)
+    delta_ii_s2_mrcc_zmq(:,1) -= delta_ij_s2_mrcc_zmq(:, i, 1) / psi_ref_coef(1,1) * psi_non_ref_coef(i, 1)
+  end do
+END_PROVIDER
 
 
  BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ]
@@ -336,8 +952,33 @@ end
   integer                        :: i, j, i_state
   
   !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
-  
-  if(mrmode == 3) then
+  if(mrmode == 4) then
+    do i = 1, N_det_ref
+      do i_state = 1, N_states
+        delta_ii(i_state,i)= delta_ii_mrcc_sto(i_state,i)
+        delta_ii_s2(i_state,i)= delta_ii_s2_mrcc_sto(i_state,i)
+      enddo
+      do j = 1, N_det_non_ref
+        do i_state = 1, N_states
+          delta_ij(i_state,j,i) = delta_ij_mrcc_sto(i_state,j,i)
+          delta_ij_s2(i_state,j,i) = delta_ij_s2_mrcc_sto(i_state,j,i)
+        enddo
+      end do
+    end do
+  else if(mrmode == 5) then
+    do i = 1, N_det_ref
+      do i_state = 1, N_states
+        delta_ii(i_state,i)= delta_ii_mrcc_zmq(i_state,i)
+        delta_ii_s2(i_state,i)= delta_ii_s2_mrcc_zmq(i_state,i)
+      enddo
+      do j = 1, N_det_non_ref
+        do i_state = 1, N_states
+          delta_ij(i_state,j,i) = delta_ij_mrcc_zmq(i_state,j,i)
+          delta_ij_s2(i_state,j,i) = delta_ij_s2_mrcc_zmq(i_state,j,i)
+        enddo
+      end do
+    end do
+  else  if(mrmode == 3) then
     do i = 1, N_det_ref
       do i_state = 1, N_states
         delta_ii(i_state,i)= delta_ii_mrcc(i_state,i)
@@ -407,6 +1048,19 @@ end
   else
     stop "invalid mrmode"
   end if
+
+  !if(mrmode == 2 .or. mrmode == 3) then
+  !   do i = 1, N_det_ref
+  !    do i_state = 1, N_states
+  !      delta_ii(i_state,i) += delta_ii_cancel(i_state,i)
+  !    enddo
+  !    do j = 1, N_det_non_ref
+  !      do i_state = 1, N_states
+  !        delta_ij(i_state,j,i) += delta_ij_cancel(i_state,j,i)
+  !      enddo
+  !    end do
+  !  end do
+  !end if
 END_PROVIDER
 
 
@@ -774,8 +1428,8 @@ end subroutine
           notf = notf+1
 
 !          call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk)
-          contrib = delta_cas(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state)
-          contrib_s2 = delta_cas_s2(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state)
+          contrib = delta_cas(II, J, i_state)* dij(J, det_cepa0_idx(k), i_state)
+          contrib_s2 = delta_cas_s2(II, J, i_state)  * dij(J, det_cepa0_idx(k), i_state)
           
           if(dabs(psi_ref_coef(J,i_state)).ge.1.d-3) then
             contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state)
diff --git a/plugins/mrcepa0/energy.irp.f b/plugins/mrcepa0/energy.irp.f
new file mode 100644
index 00000000..7eaf5562
--- /dev/null
+++ b/plugins/mrcepa0/energy.irp.f
@@ -0,0 +1,23 @@
+BEGIN_PROVIDER [ logical, initialize_mrcc_E0_denominator ]
+ implicit none
+ BEGIN_DOC
+ ! If true, initialize mrcc_E0_denominator
+ END_DOC
+ initialize_mrcc_E0_denominator = .True.
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, mrcc_E0_denominator, (N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! E0 in the denominator of the mrcc
+ END_DOC
+ if (initialize_mrcc_E0_denominator) then
+  mrcc_E0_denominator(1:N_states) = psi_energy(1:N_states)
+! mrcc_E0_denominator(1:N_states) = HF_energy - nuclear_repulsion
+! mrcc_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
+  call write_double(6,mrcc_E0_denominator(1)+nuclear_repulsion, 'mrcc Energy denominator')
+ else
+   mrcc_E0_denominator = -huge(1.d0)
+ endif
+END_PROVIDER
+
diff --git a/plugins/mrcepa0/ezfio_interface.irp.f b/plugins/mrcepa0/ezfio_interface.irp.f
new file mode 100644
index 00000000..3627abe6
--- /dev/null
+++ b/plugins/mrcepa0/ezfio_interface.irp.f
@@ -0,0 +1,80 @@
+! DO NOT MODIFY BY HAND
+! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
+! from file /home/garniron/quantum_package/src/mrcepa0/EZFIO.cfg
+
+
+BEGIN_PROVIDER [ logical, perturbative_triples  ]
+  implicit none
+  BEGIN_DOC
+! Compute perturbative contribution of the Triples
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  
+  call ezfio_has_mrcepa0_perturbative_triples(has)
+  if (has) then
+    call ezfio_get_mrcepa0_perturbative_triples(perturbative_triples)
+  else
+    print *, 'mrcepa0/perturbative_triples not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, thresh_dressed_ci  ]
+  implicit none
+  BEGIN_DOC
+! Threshold on the convergence of the dressed CI energy
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  
+  call ezfio_has_mrcepa0_thresh_dressed_ci(has)
+  if (has) then
+    call ezfio_get_mrcepa0_thresh_dressed_ci(thresh_dressed_ci)
+  else
+    print *, 'mrcepa0/thresh_dressed_ci not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_it_max_dressed_ci  ]
+  implicit none
+  BEGIN_DOC
+! Maximum number of dressed CI iterations
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  
+  call ezfio_has_mrcepa0_n_it_max_dressed_ci(has)
+  if (has) then
+    call ezfio_get_mrcepa0_n_it_max_dressed_ci(n_it_max_dressed_ci)
+  else
+    print *, 'mrcepa0/n_it_max_dressed_ci not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, lambda_type  ]
+  implicit none
+  BEGIN_DOC
+! lambda type
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  
+  call ezfio_has_mrcepa0_lambda_type(has)
+  if (has) then
+    call ezfio_get_mrcepa0_lambda_type(lambda_type)
+  else
+    print *, 'mrcepa0/lambda_type not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER
diff --git a/plugins/mrcepa0/mrcc_slave.irp.f b/plugins/mrcepa0/mrcc_slave.irp.f
new file mode 100644
index 00000000..83295985
--- /dev/null
+++ b/plugins/mrcepa0/mrcc_slave.irp.f
@@ -0,0 +1,76 @@
+program mrcc_slave
+  implicit none
+  BEGIN_DOC
+! Helper program to compute the mrcc in distributed mode.
+  END_DOC
+  print  *, "SLAVE"
+  read_wf = .False.
+  distributed_davidson = .False.
+  SOFT_TOUCH read_wf distributed_davidson
+  call provide_everything
+  call switch_qp_run_to_master
+  call run_wf
+end
+
+subroutine provide_everything
+  PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
+end
+
+subroutine run_wf
+  use f77_zmq
+  implicit none
+
+  integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
+  integer(ZMQ_PTR) :: zmq_to_qp_run_socket
+  double precision :: energy(N_states_diag)
+  character*(64) :: states(1)
+  integer :: rc, i
+  
+  call provide_everything
+  
+  zmq_context = f77_zmq_ctx_new ()
+  states(1) = 'mrcc'
+
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
+
+  do
+
+    call wait_for_states(states,zmq_state,1)
+
+    if(trim(zmq_state) == 'Stopped') then
+
+      exit
+
+    else if (trim(zmq_state) == 'mrcc') then
+
+      ! Selection
+      ! ---------
+
+      print *,  'mrcc'
+      call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
+
+      PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
+      PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
+      PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
+      PROVIDE psi_bilinear_matrix_transp_order
+
+      !$OMP PARALLEL PRIVATE(i)
+      i = omp_get_thread_num()
+      call mrcc_slave_tcp(i, energy)
+      !$OMP END PARALLEL
+      print *,  'mrcc done'
+
+    endif
+
+  end do
+end
+
+subroutine mrcc_slave_tcp(i,energy)
+  implicit none
+  double precision, intent(in) :: energy(N_states_diag)
+  integer, intent(in)            :: i
+  logical :: lstop
+  lstop = .False.
+  call run_mrcc_slave(0,i,energy,lstop)
+end
+
diff --git a/plugins/mrcepa0/mrcc_sto.irp.f b/plugins/mrcepa0/mrcc_sto.irp.f
new file mode 100644
index 00000000..5c4d318d
--- /dev/null
+++ b/plugins/mrcepa0/mrcc_sto.irp.f
@@ -0,0 +1,27 @@
+program mrsc2sub
+  implicit none
+  double precision, allocatable  :: energy(:)
+  allocate (energy(N_states))
+  
+  !!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
+  mrmode = 4
+
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call set_generators_bitmasks_as_holes_and_particles
+  if (.True.) then
+    integer :: i,j
+    do j=1,N_states
+      do i=1,N_det
+        psi_coef(i,j) = CI_eigenvectors(i,j)
+      enddo
+    enddo
+    SOFT_TOUCH psi_coef 
+  endif
+  call run(N_states,energy)
+  if(do_pt2)then
+    call run_pt2(N_states,energy)
+  endif
+  deallocate(energy)
+end
+
diff --git a/plugins/mrcepa0/mrcc_stoch.irp.f b/plugins/mrcepa0/mrcc_stoch.irp.f
new file mode 100644
index 00000000..7a11d292
--- /dev/null
+++ b/plugins/mrcepa0/mrcc_stoch.irp.f
@@ -0,0 +1,38 @@
+program mrcc_stoch
+  implicit none
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  PROVIDE mo_bielec_integrals_in_map
+  call run
+end
+
+subroutine run
+  implicit none
+  integer                        :: i,j,k
+  logical, external              :: detEq
+  
+  double precision, allocatable  :: mrcc(:)
+  integer                        :: degree
+  integer                        :: n_det_before, to_select
+  double precision               :: threshold_davidson_in
+
+  double precision               :: E_CI_before, relative_error
+
+  allocate (mrcc(N_states))
+  mrcc = 0.d0
+  !call random_seed()
+  E_CI_before = mrcc_E0_denominator(1) + nuclear_repulsion
+  threshold_selectors = 1.d0
+  threshold_generators = 1d0 
+  relative_error = 5.d-2
+  call ZMQ_mrcc(E_CI_before, mrcc, relative_error)
+  !print *,  'Final step'
+  !print *,  'N_det    = ', N_det
+  print *,  'mrcc      = ', mrcc
+  !print *,  'E        = ', E_CI_before
+  !print *,  'E+mrcc    = ', E_CI_before+mrcc
+  !print *,  '-----'
+  !call ezfio_set_full_ci_zmq_energy_mrcc(E_CI_before+mrcc(1))
+end
+
+
diff --git a/plugins/mrcepa0/mrcc_stoch_routines.irp.f b/plugins/mrcepa0/mrcc_stoch_routines.irp.f
new file mode 100644
index 00000000..0e2d9bb5
--- /dev/null
+++ b/plugins/mrcepa0/mrcc_stoch_routines.irp.f
@@ -0,0 +1,604 @@
+BEGIN_PROVIDER [ integer, fragment_first ]
+  implicit none
+  fragment_first = first_det_of_teeth(1)
+END_PROVIDER
+
+
+subroutine ZMQ_mrcc(E, mrcc, delta, delta_s2, relative_error)
+  use dress_types
+  use f77_zmq
+  
+  implicit none
+  
+  character(len=64000)           :: task
+  integer(ZMQ_PTR)               :: zmq_to_qp_run_socket, zmq_to_qp_run_socket2
+  integer, external              :: omp_get_thread_num
+  double precision, intent(in)   :: relative_error, E
+  double precision, intent(out)  :: mrcc(N_states)
+  double precision, intent(out)  :: delta(N_states, N_det_non_ref)
+  double precision, intent(out)  :: delta_s2(N_states, N_det_non_ref)
+  
+  
+  integer                        :: i, j, k, Ncp
+  
+  double precision, external     :: omp_get_wtime
+  double precision               :: time
+
+
+  
+  provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral mrcc_weight psi_selectors
+  
+  
+  
+  print *, '========== ================= ================= ================='
+  print *, ' Samples        Energy         Stat. Error         Seconds      '
+  print *, '========== ================= ================= ================='
+  
+  call new_parallel_job(zmq_to_qp_run_socket,'mrcc')
+  call zmq_put_psi(zmq_to_qp_run_socket,1,mrcc_e0_denominator,size(mrcc_e0_denominator))
+  
+!   call get_carlo_workbatch(Ncp, tbc, cp, cp_at, cp_N)
+
+  do i=1,comb_teeth
+    print *, "TOOTH", first_det_of_teeth(i+1) - first_det_of_teeth(i)
+  end do
+  !stop
+
+  integer(ZMQ_PTR), external     :: new_zmq_to_qp_run_socket
+  integer                        :: ipos
+  ipos=1
+  do i=1,N_mrcc_jobs
+    if(mrcc_jobs(i) > fragment_first) then
+      write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') 0, mrcc_jobs(i)
+      ipos += 20
+      if (ipos > 63980) then
+        call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos)))
+        ipos=1
+      endif
+    else
+      do j=1,fragment_count
+        write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') j, mrcc_jobs(i)
+        ipos += 20
+        if (ipos > 63980) then
+          call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos)))
+          ipos=1
+        endif
+      end do
+    end if
+  end do
+  if (ipos > 1) then
+    call add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos)))
+  endif
+  
+  call zmq_set_running(zmq_to_qp_run_socket)
+  
+  !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc+1)              &
+      !$OMP  PRIVATE(i)
+  i = omp_get_thread_num()
+  if (i==0) then
+    call mrcc_collector(E, relative_error, delta, delta_s2, mrcc)
+  else
+    call mrcc_slave_inproc(i)
+  endif
+  !$OMP END PARALLEL
+  call end_parallel_job(zmq_to_qp_run_socket, 'mrcc')
+  
+  print *, '========== ================= ================= ================='
+end subroutine
+
+
+subroutine mrcc_slave_inproc(i)
+  implicit none
+  integer, intent(in)            :: i
+
+  call run_mrcc_slave(1,i,mrcc_e0_denominator)
+end
+
+
+subroutine mrcc_collector(E, relative_error, delta, delta_s2, mrcc)
+  use dress_types
+  use f77_zmq
+  use bitmasks
+  implicit none
+
+  
+  double precision, intent(in) :: relative_error, E
+  double precision, intent(out)  :: mrcc(N_states)
+  double precision, allocatable  :: cp(:,:,:,:)
+
+  double precision, intent(out)  :: delta(N_states, N_det_non_ref)
+  double precision, intent(out)  :: delta_s2(N_states, N_det_non_ref)
+  double precision, allocatable  :: delta_loc(:,:,:), delta_det(:,:,:,:)
+  double precision, allocatable  :: mrcc_detail(:,:)
+  double precision               :: mrcc_mwen(N_states)
+  integer(ZMQ_PTR),external      :: new_zmq_to_qp_run_socket
+  integer(ZMQ_PTR)               :: zmq_to_qp_run_socket
+
+  integer(ZMQ_PTR), external     :: new_zmq_pull_socket
+  integer(ZMQ_PTR)               :: zmq_socket_pull
+
+  integer :: more
+  integer :: i, j, k, i_state, N, ntask
+  integer, allocatable :: task_id(:)
+  integer :: Nindex
+  integer, allocatable :: ind(:)
+  double precision, save :: time0 = -1.d0
+  double precision :: time, timeLast, old_tooth
+  double precision, external :: omp_get_wtime
+  integer :: cur_cp, old_cur_cp
+  integer, allocatable :: parts_to_get(:)
+  logical, allocatable :: actually_computed(:)
+  integer :: total_computed
+  
+  allocate(delta_det(N_states, N_det_non_ref, 0:comb_teeth+1, 2))
+  allocate(cp(N_states, N_det_non_ref, N_cp, 2), mrcc_detail(N_states, N_det_generators))
+  allocate(delta_loc(N_states, N_det_non_ref, 2))
+  mrcc_detail = 0d0
+  delta_det = 0d0
+  !mrcc_detail = mrcc_detail / 0d0
+  cp = 0d0
+  total_computed = 0
+  character*(512) :: task
+  
+  allocate(actually_computed(N_det_generators), parts_to_get(N_det_generators))
+    
+  mrcc_mwen =0.d0
+  
+  parts_to_get(:) = 1
+  if(fragment_first > 0) then
+    do i=1,fragment_first
+      parts_to_get(i) = fragment_count
+    enddo
+  endif
+
+  actually_computed = .false.
+
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
+  zmq_socket_pull = new_zmq_pull_socket()
+  allocate(task_id(N_det_generators), ind(1))
+  more = 1
+  if (time0 < 0.d0) then
+      call wall_time(time0)
+  endif
+  timeLast = time0
+  cur_cp = 0
+  old_cur_cp = 0
+  pullLoop : do while (more == 1)
+    call pull_mrcc_results(zmq_socket_pull, Nindex, ind, mrcc_mwen, delta_loc, task_id, ntask)
+
+    if(Nindex /= 1) stop "tried pull multiple Nindex"
+
+    do i=1,Nindex
+      mrcc_detail(:, ind(i)) += mrcc_mwen(:)
+      do j=1,N_cp !! optimizable
+        if(cps(ind(i), j) > 0d0) then
+          if(tooth_of_det(ind(i)) < cp_first_tooth(j)) stop "coef on supposedely deterministic det"
+          double precision :: fac
+          integer :: toothMwen
+          logical :: fracted
+          fac = cps(ind(i), j) / cps_N(j) * mrcc_weight_inv(ind(i)) * comb_step
+          do k=1,N_det_non_ref
+          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
+        end if
+      end do
+      toothMwen = tooth_of_det(ind(i))
+      fracted = (toothMwen /= 0)
+      if(fracted) fracted = (ind(i) == 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))
+      else
+        delta_det(:,:,toothMwen, 1) += delta_loc(:,:,1)
+        delta_det(:,:,toothMwen, 2) += delta_loc(:,:,2)
+      end if
+
+      parts_to_get(ind(i)) -= 1
+      if(parts_to_get(ind(i)) == 0) then
+        actually_computed(ind(i)) = .true.
+        !print *, "CONTRIB", ind(i), psi_non_ref_coef(ind(i),1), mrcc_detail(1, ind(i))
+        total_computed += 1
+      end if
+    end do
+
+    do i=1, ntask
+      call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more)
+    end do
+     
+    time = omp_get_wtime()
+    
+    
+
+    if(time - timeLast > 1d0 .or. more /= 1) then
+      timeLast = time
+      cur_cp = N_cp
+      if(.not. actually_computed(mrcc_jobs(1))) cycle pullLoop
+
+      do i=2,N_det_generators
+        if(.not. actually_computed(mrcc_jobs(i))) then
+          print *, "first not comp", i
+          cur_cp = done_cp_at(i-1)
+          exit
+        end if
+      end do
+      if(cur_cp == 0) cycle pullLoop
+      
+      !!!!!!!!!!!!
+      double precision :: su, su2, eqt, avg, E0, val
+      su = 0d0
+      su2 = 0d0
+      
+      if(N_states > 1) stop "mrcc_stoch : N_states == 1"
+      do i=1, int(cps_N(cur_cp))
+        !if(.not. actually_computed(i)) stop "not computed"
+        !call get_comb_val(comb(i), mrcc_detail, cp_first_tooth(cur_cp), val)
+        call get_comb_val(comb(i), mrcc_detail, cur_cp, val)
+        !val = mrcc_detail(1, i) * mrcc_weight_inv(i) * comb_step
+        su += val ! cps(i, cur_cp) * val
+        su2 += val**2 ! cps(i, cur_cp) * val**2
+      end do
+      avg = su / cps_N(cur_cp)
+      eqt = dsqrt( ((su2 / cps_N(cur_cp)) - avg**2) / cps_N(cur_cp) )
+      E0 = sum(mrcc_detail(1, :first_det_of_teeth(cp_first_tooth(cur_cp))-1))
+      if(cp_first_tooth(cur_cp) <= comb_teeth) then
+        E0 = E0 + mrcc_detail(1, 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
+        ! 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
+        call zmq_abort(zmq_to_qp_run_socket)
+      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, ''
+        endif
+      endif
+    end if
+  end do pullLoop
+ 
+  if(total_computed == N_det_generators) then
+    print *, "TOTALLY COMPUTED"
+    delta = 0d0
+    delta_s2 = 0d0
+    do i=comb_teeth+1,0,-1
+      delta += delta_det(:,:,i,1)
+      delta_s2 += delta_det(:,:,i,2)
+    end do
+  else
+
+
+  delta = cp(:,:,cur_cp,1)
+  delta_s2 = cp(:,:,cur_cp,2)
+  
+  do i=cp_first_tooth(cur_cp)-1,0,-1
+    delta += delta_det(:,:,i,1)
+    delta_s2 += delta_det(:,:,i,2)
+  end do
+
+  end if
+
+  mrcc(1) = E
+  
+  call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
+  call end_zmq_pull_socket(zmq_socket_pull)
+end subroutine
+
+
+integer function mrcc_find(v, w, sze, imin, imax)
+  implicit none
+  integer, intent(in) :: sze, imin, imax
+  double precision, intent(in) :: v, w(sze)
+  integer :: i,l,h
+  integer, parameter :: block=64
+
+  l = imin
+  h = imax-1
+
+  do while(h-l >= block)
+    i = ishft(h+l,-1)
+    if(w(i+1) > v) then
+      h = i-1
+    else
+      l = i+1
+    end if
+  end do
+  !DIR$ LOOP COUNT (64)
+  do mrcc_find=l,h
+    if(w(mrcc_find) >= v) then
+      exit
+    end if
+  end do
+end function
+
+
+ BEGIN_PROVIDER [ integer, gen_per_cp ]
+&BEGIN_PROVIDER [ integer, comb_teeth ]
+&BEGIN_PROVIDER [ integer, N_cps_max ]
+  implicit none
+  comb_teeth = 16
+  N_cps_max = 32
+  !comb_per_cp = 64
+  gen_per_cp = (N_det_generators / N_cps_max) + 1
+  N_cps_max += 1
+  !N_cps_max = N_det_generators / comb_per_cp + 1
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer, N_cp ]
+&BEGIN_PROVIDER [ double precision, cps_N, (N_cps_max) ]
+&BEGIN_PROVIDER [ integer, cp_first_tooth, (N_cps_max) ]
+&BEGIN_PROVIDER [ integer, done_cp_at, (N_det_generators) ]
+&BEGIN_PROVIDER [ double precision, cps, (N_det_generators, N_cps_max) ]
+&BEGIN_PROVIDER [ integer, N_mrcc_jobs ]
+&BEGIN_PROVIDER [ integer, mrcc_jobs, (N_det_generators) ]
+&BEGIN_PROVIDER [ double precision, comb, (N_det_generators) ]
+! subroutine get_carlo_workbatch(Ncp, tbc, cps, done_cp_at)
+  implicit none
+  logical, allocatable         :: computed(:)
+  integer                        :: i, j, last_full, dets(comb_teeth)
+  integer                        :: k, l, cur_cp, under_det(comb_teeth+1)
+  integer, allocatable :: iorder(:), first_cp(:)
+
+  allocate(iorder(N_det_generators), first_cp(N_cps_max+1))
+  allocate(computed(N_det_generators))
+  first_cp = 1
+  cps = 0d0
+  cur_cp = 1
+  done_cp_at = 0
+  
+  computed = .false.
+  
+  N_mrcc_jobs = first_det_of_comb - 1
+  do i=1, N_mrcc_jobs
+    mrcc_jobs(i) = i
+    computed(i) = .true.
+  end do
+  
+  l=first_det_of_comb
+  call RANDOM_NUMBER(comb)
+  do i=1,N_det_generators
+    comb(i) = comb(i) * comb_step
+    !DIR$ FORCEINLINE
+    call add_comb(comb(i), computed, cps(1, cur_cp), N_mrcc_jobs, mrcc_jobs)
+    
+    if(N_mrcc_jobs / gen_per_cp > (cur_cp-1) .or. N_mrcc_jobs == N_det_generators) then
+    !if(mod(i, comb_per_cp) == 0 .or. N_mrcc_jobs == N_det_generators) then
+      first_cp(cur_cp+1) = N_mrcc_jobs
+      done_cp_at(N_mrcc_jobs) = cur_cp
+      cps_N(cur_cp) = dfloat(i)
+      if(N_mrcc_jobs /= N_det_generators) then
+        cps(:, cur_cp+1) = cps(:,  cur_cp)
+        cur_cp += 1
+      end if
+      !cps(:, cur_cp) = cps(:, cur_cp) / dfloat(i)
+      
+      if (N_mrcc_jobs == N_det_generators) exit
+    end if
+    do while (computed(l))
+      l=l+1
+    enddo
+    k=N_mrcc_jobs+1
+    mrcc_jobs(k) = l
+    computed(l) = .True.
+    N_mrcc_jobs = k
+  enddo
+  N_cp = cur_cp
+  if(N_mrcc_jobs /= N_det_generators .or. N_cp > N_cps_max) then
+    print *, N_mrcc_jobs, N_det_generators, N_cp, N_cps_max
+    stop "carlo workcarlo_workbatch"
+  end if
+
+  cur_cp = 0
+  do i=1,N_mrcc_jobs
+    if(done_cp_at(i) /= 0) cur_cp = done_cp_at(i)
+    done_cp_at(i) = cur_cp
+  end do
+  
+
+  under_det = 0
+  cp_first_tooth = 0
+  do i=1,N_mrcc_jobs
+    do j=comb_teeth+1,1,-1
+      if(mrcc_jobs(i) <= first_det_of_teeth(j)) then
+        under_det(j) = under_det(j) + 1
+        if(under_det(j) == first_det_of_teeth(j))then
+          do l=done_cp_at(i)+1, N_cp
+            cps(:first_det_of_teeth(j)-1, l) = 0d0
+            cp_first_tooth(l) = j
+          end do
+          cps(first_det_of_teeth(j), done_cp_at(i)+1) = &
+              cps(first_det_of_teeth(j), done_cp_at(i)+1) * fractage(j)
+        end if
+      else
+        exit
+      end if
+    end do
+  end do
+  cps(:, N_cp) = 0d0
+  cp_first_tooth(N_cp) = comb_teeth+1
+
+  iorder = -1132154665
+  do i=1,N_cp-1
+    call isort(mrcc_jobs(first_cp(i)+1:first_cp(i+1)),iorder,first_cp(i+1)-first_cp(i))
+  end do
+! end subroutine
+END_PROVIDER
+
+
+subroutine get_comb_val(stato, detail, cur_cp, val)
+  implicit none
+  integer, intent(in)           :: cur_cp
+  integer                       :: first
+  double precision, intent(in)  :: stato, detail(N_states, N_det_generators)
+  double precision, intent(out) :: val
+  double precision :: curs
+  integer :: j, k
+  integer, external :: mrcc_find
+
+  curs = 1d0 - stato
+  val = 0d0
+  first = cp_first_tooth(cur_cp) 
+
+  do j = comb_teeth, first, -1
+    !DIR$ FORCEINLINE
+    k = mrcc_find(curs, mrcc_cweight,size(mrcc_cweight), first_det_of_teeth(j), first_det_of_teeth(j+1))
+    !if(k < first_det_of_teeth(first)) exit
+    if(k == first_det_of_teeth(first)) then
+     val += detail(1, k) * mrcc_weight_inv(k) * comb_step * fractage(first)
+    else
+     val += detail(1, k) * mrcc_weight_inv(k) * comb_step
+    end if
+
+    curs -= comb_step
+  end do
+end subroutine
+
+
+subroutine get_comb(stato, dets)
+  implicit none
+  double precision, intent(in) :: stato
+  integer, intent(out) :: dets(comb_teeth)
+  double precision :: curs
+  integer :: j
+  integer, external :: mrcc_find
+
+  curs = 1d0 - stato
+  do j = comb_teeth, 1, -1
+    !DIR$ FORCEINLINE
+    dets(j) = mrcc_find(curs, mrcc_cweight,size(mrcc_cweight), first_det_of_teeth(j), first_det_of_teeth(j+1))
+    curs -= comb_step
+  end do
+end subroutine
+
+
+subroutine add_comb(com, computed, cp, N, tbc)
+  implicit none
+  double precision, intent(in) :: com
+  integer, intent(inout) :: N
+  double precision, intent(inout) :: cp(N_det_non_ref)
+  logical, intent(inout) :: computed(N_det_generators)
+  integer, intent(inout) :: tbc(N_det_generators)
+  integer :: i, k, l, dets(comb_teeth)
+  
+  !DIR$ FORCEINLINE
+  call get_comb(com, dets)
+  
+  k=N+1
+  do i = 1, comb_teeth
+    l = dets(i)
+    cp(l) += 1d0 ! mrcc_weight_inv(l) * comb_step
+    if(.not.(computed(l))) then
+      tbc(k) = l
+      k = k+1
+      computed(l) = .true.
+    end if
+  end do
+  N = k-1
+end subroutine
+
+
+
+ BEGIN_PROVIDER [ double precision, mrcc_weight, (N_det_generators) ]
+&BEGIN_PROVIDER [ double precision, mrcc_weight_inv, (N_det_generators) ]
+&BEGIN_PROVIDER [ double precision, mrcc_cweight, (N_det_generators) ]
+&BEGIN_PROVIDER [ double precision, mrcc_cweight_cache, (N_det_generators) ]
+&BEGIN_PROVIDER [ double precision, fractage, (comb_teeth) ]
+&BEGIN_PROVIDER [ double precision, comb_step ]
+&BEGIN_PROVIDER [ integer, first_det_of_teeth, (comb_teeth+1) ]
+&BEGIN_PROVIDER [ integer, first_det_of_comb ]
+&BEGIN_PROVIDER [ integer, tooth_of_det, (N_det_generators) ]
+  implicit none
+  integer :: i
+  double precision :: norm_left, stato
+  integer, external :: mrcc_find  
+
+  mrcc_weight(1) = psi_coef_generators(1,1)**2
+  mrcc_cweight(1) = psi_coef_generators(1,1)**2
+  
+  do i=1,N_det_generators
+    mrcc_weight(i) = psi_coef_generators(i,1)**2
+  enddo
+
+  ! Important to loop backwards for numerical precision
+  mrcc_cweight(N_det_generators) = mrcc_weight(N_det_generators)
+  do i=N_det_generators-1,1,-1
+    mrcc_cweight(i) = mrcc_weight(i) + mrcc_cweight(i+1) 
+  end do
+  
+  do i=1,N_det_generators
+    mrcc_weight(i)  = mrcc_weight(i) / mrcc_cweight(1)
+    mrcc_cweight(i) = mrcc_cweight(i) / mrcc_cweight(1)
+  enddo
+
+  do i=1,N_det_generators-1
+    mrcc_cweight(i) = 1.d0 - mrcc_cweight(i+1) 
+  end do
+  mrcc_cweight(N_det_generators) = 1.d0
+  
+  norm_left = 1d0
+  
+  comb_step = 1d0/dfloat(comb_teeth)
+  first_det_of_comb = 1
+  do i=1,N_det_generators
+    if(mrcc_weight(i)/norm_left < .25d0*comb_step) then
+      first_det_of_comb = i
+      exit
+    end if
+    norm_left -= mrcc_weight(i)
+  end do
+  first_det_of_comb = max(2,first_det_of_comb)
+  call write_int(6, first_det_of_comb-1, 'Size of deterministic set')
+  
+
+  comb_step =  (1d0 - mrcc_cweight(first_det_of_comb-1)) * comb_step
+  
+  stato = 1d0 - comb_step
+  iloc = N_det_generators
+  do i=comb_teeth, 1, -1
+    integer :: iloc
+    iloc = mrcc_find(stato, mrcc_cweight, N_det_generators, 1, iloc)
+    first_det_of_teeth(i) = iloc
+    fractage(i) = (mrcc_cweight(iloc) - stato) / mrcc_weight(iloc)
+    stato -= comb_step
+  end do
+  first_det_of_teeth(comb_teeth+1) = N_det_generators + 1
+  first_det_of_teeth(1) = first_det_of_comb
+  
+
+  if(first_det_of_teeth(1) /= first_det_of_comb) then
+     print *, 'Error in ', irp_here
+     stop "comb provider"
+  endif
+  
+  do i=1,N_det_generators
+    mrcc_weight_inv(i) = 1.d0/mrcc_weight(i)
+  enddo
+
+  tooth_of_det(:first_det_of_teeth(1)-1) = 0
+  do i=1,comb_teeth
+    tooth_of_det(first_det_of_teeth(i):first_det_of_teeth(i+1)-1) = i
+  end do
+
+  !double precision :: cur
+  !fractage = 1d0
+  !do i=1,comb_teeth-1
+  !  cur = 1d0 - dfloat(i)*comb_step
+    
+  !end do
+END_PROVIDER
+
+
+
+
+
+
+
diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f
index 53fa74a1..4bb225a3 100644
--- a/plugins/mrcepa0/mrcepa0_general.irp.f
+++ b/plugins/mrcepa0/mrcepa0_general.irp.f
@@ -41,12 +41,12 @@ subroutine run(N_st,energy)
       print *,  'Iteration', iteration, '/', n_it_mrcc_max
       print *,  '===============================================' 
       print *,  ''
-      E_old = sum(ci_energy_dressed(1:N_states))
+      E_old = mrcc_e0_denominator(1) !sum(ci_energy_dressed(1:N_states))
       do i=1,N_st
         call write_double(6,ci_energy_dressed(i),"Energy")
       enddo
       call diagonalize_ci_dressed(lambda)
-      E_new = sum(ci_energy_dressed(1:N_states))
+      E_new = mrcc_e0_denominator(1) !sum(ci_energy_dressed(1:N_states))
 !  if (.true.) then
 !    provide delta_ij_mrcc_pouet
 !  endif
diff --git a/plugins/mrcepa0/run_mrcc_slave.irp.f b/plugins/mrcepa0/run_mrcc_slave.irp.f
new file mode 100644
index 00000000..b6088424
--- /dev/null
+++ b/plugins/mrcepa0/run_mrcc_slave.irp.f
@@ -0,0 +1,244 @@
+
+BEGIN_PROVIDER [ integer, fragment_count ]
+  implicit none
+  BEGIN_DOC
+  ! Number of fragments for the deterministic part
+  END_DOC
+  fragment_count = 1 !! (elec_alpha_num-n_core_orb)**2
+END_PROVIDER
+
+
+subroutine run_mrcc_slave(thread,iproc,energy)
+  use dress_types
+  use f77_zmq
+  implicit none
+
+  double precision, intent(in)    :: energy(N_states_diag)
+  integer,  intent(in)            :: thread, iproc
+  integer                        :: rc, i
+
+  integer                        :: worker_id, task_id(1), ctask, ltask
+  character*(512)                :: task
+
+  integer(ZMQ_PTR),external      :: new_zmq_to_qp_run_socket
+  integer(ZMQ_PTR)               :: zmq_to_qp_run_socket
+
+  integer(ZMQ_PTR), external     :: new_zmq_push_socket
+  integer(ZMQ_PTR)               :: zmq_socket_push
+
+  logical :: done
+
+  double precision,allocatable :: mrcc_detail(:,:)
+  integer :: ind(1)
+  integer :: Nindex
+  
+  integer(bit_kind),allocatable :: abuf(:,:,:)
+  integer(bit_kind) :: mask(N_int,2), omask(N_int,2)
+ 
+  double precision,allocatable :: delta_ij_loc(:,:,:) 
+  double precision,allocatable :: delta_ii_loc(:,:) 
+  !double precision,allocatable :: delta_ij_s2_loc(:,:,:)
+  !double precision,allocatable :: delta_ii_s2_loc(:,:)
+  integer :: h,p,n
+  logical :: ok
+  double precision :: contrib(N_states)
+
+  allocate(delta_ij_loc(N_states,N_det_non_ref,2) &
+  ,delta_ii_loc(N_states,2))! &
+  !,delta_ij_s2_loc(N_states,N_det_non_ref,N_det_ref) &
+  !,delta_ii_s2_loc(N_states, N_det_ref))
+
+
+        allocate(abuf(N_int, 2, N_det_non_ref))
+
+  allocate(mrcc_detail(N_states, N_det_generators))
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
+  zmq_socket_push      = new_zmq_push_socket(thread)
+  call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
+  if(worker_id == -1) then
+    print *, "WORKER -1"
+    call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
+    call end_zmq_push_socket(zmq_socket_push,thread)
+    return
+  end if
+  !buf%N = 0
+  ctask = 1
+  Nindex=1
+  mrcc_detail = 0d0
+  do
+    call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id(ctask), task)
+
+    done = task_id(ctask) == 0
+    if (done) then
+      ctask = ctask - 1
+    else
+      integer :: i_generator, i_i_generator, subset
+      read (task,*) subset, ind
+      
+!      if(buf%N == 0) then
+!        ! Only first time 
+!        call create_selection_buffer(1, 2, buf)
+!      end if
+      do i_i_generator=1, Nindex
+        contrib = 0d0
+        i_generator = ind(i_i_generator)
+        delta_ij_loc = 0d0
+        delta_ii_loc = 0d0
+        !delta_ij_s2_loc = 0d0
+        !delta_ii_s2_loc = 0d0
+        !call select_connected(i_generator,energy,mrcc_detail(1, i_i_generator),buf,subset)
+
+!!!!!!!!!!!!!!!!!!!!!!
+      !if(mod(gen, 1000) == 0) print *, "mrcc ", gen, "/", N_det_generators
+        do h=1, hh_shortcut(0)
+          call apply_hole_local(psi_det_generators(1,1,i_generator), hh_exists(1, h), mask, ok, N_int)
+          if(.not. ok) cycle
+          omask = 0_bit_kind
+          if(hh_exists(1, h) /= 0) omask = mask
+          n = 1
+          do p=hh_shortcut(h), hh_shortcut(h+1)-1
+            call apply_particle_local(mask, pp_exists(1, p), abuf(1,1,n), ok, N_int)
+            if(ok) n = n + 1
+            if(n > N_det_non_ref) stop "Buffer too small in MRCC..."
+          end do
+          n = n - 1
+
+          if(n /= 0) then
+            call mrcc_part_dress_1c(delta_ij_loc(1,1,1), delta_ii_loc(1,1), delta_ij_loc(1,1,2), delta_ii_loc(1,2), &
+                  i_generator,n,abuf,N_int,omask,contrib)
+          endif
+        end do
+!!!!!!!!!!!!!!!!!!!!!!
+        !print *, "contrib", i_generator, contrib
+        mrcc_detail(:, i_i_generator) = contrib
+        if(Nindex /= 1) stop "run_mrcc_slave multiple dress not supported"
+      enddo
+    endif
+
+    if(done .or. (ctask == size(task_id)) ) then
+!      if(buf%N == 0 .and. ctask > 0) stop "uninitialized selection_buffer"
+      do i=1, ctask
+         call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i))
+      end do
+      if(ctask > 0) then
+        call push_mrcc_results(zmq_socket_push, Nindex, ind, mrcc_detail, delta_ij_loc(1,1,1), task_id(1), ctask)
+        mrcc_detail(:,:Nindex) = 0d0
+!        buf%cur = 0
+      end if
+      ctask = 0
+    end if
+
+    if(done) exit
+    ctask = ctask + 1
+  end do
+  call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
+  call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
+  call end_zmq_push_socket(zmq_socket_push,thread)
+!  call delete_selection_buffer(buf)
+end subroutine
+
+
+subroutine push_mrcc_results(zmq_socket_push, N, ind, mrcc_detail, delta_loc, task_id, ntask)
+  use f77_zmq
+  implicit none
+
+  integer(ZMQ_PTR), intent(in)   :: zmq_socket_push
+  double precision, intent(in)   :: mrcc_detail(N_states, N_det_generators)
+  double precision, intent(in)   :: delta_loc(N_states, N_det_non_ref, 2)
+  integer, intent(in) :: ntask, N, ind(*), task_id(*)
+  integer :: rc, i
+  
+  if(N/=1) stop "mrcc_stoch, tried to push N>1"
+
+  rc = f77_zmq_send( zmq_socket_push, N, 4, ZMQ_SNDMORE)
+  if(rc /= 4) stop "push"
+
+  rc = f77_zmq_send( zmq_socket_push, ind, 4*N, ZMQ_SNDMORE)
+  if(rc /= 4*N) stop "push"
+
+
+  rc = f77_zmq_send( zmq_socket_push, mrcc_detail, 8*N_states*N, ZMQ_SNDMORE)
+  if(rc /= 8*N_states*N) stop "push"
+
+  rc = f77_zmq_send( zmq_socket_push, delta_loc(1,1,1), 8*N_states*N_det_non_ref*N, ZMQ_SNDMORE)
+  if(rc /= 8*N_states*N_det_non_ref*N) stop "push"
+  
+  rc = f77_zmq_send( zmq_socket_push, delta_loc(1,1,2), 8*N_states*N_det_non_ref*N, ZMQ_SNDMORE)
+  if(rc /= 8*N_states*N_det_non_ref*N) stop "push"
+
+
+  rc = f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE)
+  if(rc /= 4) stop "push"
+
+  rc = f77_zmq_send( zmq_socket_push, task_id, ntask*4, 0)
+  if(rc /= 4*ntask) stop "push"
+
+! Activate is zmq_socket_push is a REQ
+IRP_IF ZMQ_PUSH
+IRP_ELSE
+  character*(2) :: ok
+  rc = f77_zmq_recv( zmq_socket_push, ok, 2, 0)
+IRP_ENDIF
+
+end subroutine
+
+
+subroutine pull_mrcc_results(zmq_socket_pull, N, ind, mrcc_detail, delta_loc, task_id, ntask)
+  use dress_types
+  use f77_zmq
+  implicit none
+  integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull
+  double precision, intent(inout) :: mrcc_detail(N_states)
+  double precision, intent(inout) :: delta_loc(N_states, N_det_non_ref, 2)
+  double precision, allocatable :: mrcc_dress_mwen(:,:)
+  integer, intent(out) :: ind(*)
+  integer, intent(out) :: N, ntask, task_id(*)
+  integer :: rc, rn, i
+
+  allocate(mrcc_dress_mwen(N_states, N_det_non_ref))
+  
+  rc = f77_zmq_recv( zmq_socket_pull, N, 4, 0)
+  if(rc /= 4) stop "pull"
+  if(N /= 1) stop "mrcc : pull with N>1 not supported"
+
+  rc = f77_zmq_recv( zmq_socket_pull, ind, 4, 0)
+  if(rc /= 4*N) stop "pull"
+  
+  rc = f77_zmq_recv( zmq_socket_pull, mrcc_detail, N_states*8*N, 0)
+  if(rc /= 8*N_states*N) stop "pull"
+  
+  !do i=1,N
+    rc = f77_zmq_recv( zmq_socket_pull, delta_loc(1,1,1), N_states*8*N_det_non_ref, 0)
+    if(rc /= 8*N_states*N_det_non_ref) stop "pull"
+    !call add_to_dress_buffer(mrcc_dress(1), ind(i), mrcc_dress_mwen)
+  !end do
+  
+  !do i=1,N
+    rc = f77_zmq_recv( zmq_socket_pull, delta_loc(1,1,2), N_states*8*N_det_non_ref, 0)
+    if(rc /= 8*N_states*N_det_non_ref) stop "pull"
+    !call add_to_dress_buffer(mrcc_dress(2), ind(i), mrcc_dress_mwen)
+  !end do
+
+  rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)
+  if(rc /= 4) stop "pull"
+
+  rc = f77_zmq_recv( zmq_socket_pull, task_id, ntask*4, 0)
+  if(rc /= 4*ntask) stop "pull"
+
+! Activate is zmq_socket_pull is a REP
+IRP_IF ZMQ_PUSH
+IRP_ELSE
+  rc = f77_zmq_send( zmq_socket_pull, 'ok', 2, 0)
+IRP_ENDIF
+
+end subroutine
+ 
+ 
+BEGIN_PROVIDER [ double precision, mrcc_workload, (N_det_generators) ]
+  integer :: i
+  do i=1,N_det_generators
+    mrcc_workload(i) = dfloat(N_det_generators - i + 1)**2
+  end do
+  mrcc_workload = mrcc_workload / sum(mrcc_workload)
+END_PROVIDER
+