diff --git a/plugins/Psiref_threshold/psi_ref.irp.f b/plugins/Psiref_threshold/psi_ref.irp.f
index 5e722822..ee69ef5c 100644
--- a/plugins/Psiref_threshold/psi_ref.irp.f
+++ b/plugins/Psiref_threshold/psi_ref.irp.f
@@ -6,19 +6,22 @@ use bitmasks
 &BEGIN_PROVIDER [ integer, N_det_ref ]
   implicit none
   BEGIN_DOC
-  ! Reference wave function, defined as determinants with coefficients > 0.05
+  ! Reference wave function, defined as determinants with amplitudes > 0.05
   ! idx_ref gives the indice of the ref determinant in psi_det.
   END_DOC
   integer                        :: i, k, l
   logical                        :: good
-  double precision, parameter :: threshold=0.05d0
+  double precision, parameter    :: threshold=0.05d0
+  double precision               :: t(N_states)
   N_det_ref = 0
-  t = threshold * abs_psi_coef_max
+  do l = 1, N_states
+    t(l) = threshold * abs_psi_coef_max(l)
+  enddo
   do i=1,N_det
     good = .False.
-    do l = 1, N_states
+    do l=1, N_states
      psi_ref_coef(i,l) = 0.d0
-     good = good.or.(dabs(psi_coef(i,l)) > t)
+     good = good.or.(dabs(psi_coef(i,l)) > t(l))
     enddo
     if (good) then
       N_det_ref = N_det_ref+1
diff --git a/plugins/mrcc_selected/dressing.irp.f b/plugins/mrcc_selected/dressing.irp.f
new file mode 100644
index 00000000..3646b0b2
--- /dev/null
+++ b/plugins/mrcc_selected/dressing.irp.f
@@ -0,0 +1,1008 @@
+use bitmasks
+
+
+
+ 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) ]
+  use bitmasks
+  implicit none
+  integer :: gen, h, p, n, t, i, 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
+  
+  delta_ij_mrcc = 0d0
+  delta_ii_mrcc = 0d0
+  print *, "Dij", dij(1,1,1)
+  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(N_det_non_ref, N_det_ref, delta_ii_mrcc, delta_ij_mrcc) &
+  !$OMP private(h, n, mask, omask, buf, ok, iproc)
+  do gen= 1, N_det_generators
+    allocate(buf(N_int, 2, N_det_non_ref))
+    iproc = omp_get_thread_num() + 1
+    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,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 "MRCC..."
+      end do
+      n = n - 1
+
+      if(n /= 0) call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc,gen,n,buf,N_int,omask)
+
+    end do
+    deallocate(buf)
+  end do
+  !$OMP END PARALLEL DO
+END_PROVIDER
+
+
+! subroutine blit(b1, b2)
+!   double precision :: b1(N_states,N_det_non_ref,N_det_ref), b2(N_states,N_det_non_ref,N_det_ref)
+!   b1 = b1 + b2
+! end subroutine
+
+
+subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffer,Nint,key_mask)
+ use bitmasks
+ implicit none
+
+  integer, intent(in)            :: i_generator,n_selected, Nint
+  double precision, intent(inout) :: delta_ij_(N_states,N_det_non_ref,N_det_ref)
+  double precision, intent(inout) :: delta_ii_(N_states,N_det_ref)
+
+  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 ,degree
+
+  double precision               :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states)
+  double precision, allocatable  :: dIa_hla(:,:)
+  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(:)
+  
+  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
+  
+
+  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))
+  allocate(idx_alpha(0:psi_det_size), degree_alpha(psi_det_size))
+  !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint)
+  call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
+  
+!   if(fullMatch) then
+!     return
+!   end if
+  
+  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))
+  
+  ! |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(key_mask(1,1) /= 0) then !!!!!!!!!!! PAS GENERAL !!!!!!!!!
+      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))
+    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),degree,Nint)
+      if (degree > 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)
+        ! Loop if lambda == 0
+        logical                        :: loop
+!         loop = .True.
+!         do i_state=1,N_states
+!           if (lambda_mrcc(i_state,idx_alpha(k_sd)) /= 0.d0) then
+!             loop = .False.
+!             exit
+!           endif
+!         enddo
+!         if (loop) then
+!           cycle
+!         endif
+        
+        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>
+        ! <I|H|k>
+        !hIk = hij_mrcc(idx_alpha(k_sd),i_I)
+        !         call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk)
+        
+        do i_state=1,N_states
+          dIK(i_state) = dij(i_I, idx_alpha(k_sd), i_state)
+          !dIk(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(k_sd)), N_int) !!hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
+          !dIk(i_state) = psi_non_ref_coef(idx_alpha(k_sd), i_state) / psi_ref_coef(i_I, i_state)
+        enddo
+        
+        
+        ! |l> = Exc(k -> alpha) |I>
+        call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
+        call decode_exc(exc,degree,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)
+        if(.not. ok) cycle
+        
+        ! <I| \l/ |alpha>
+        do i_state=1,N_states
+          dka(i_state) = 0.d0
+        enddo
+        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
+            
+!             loop = .True.
+!             do i_state=1,N_states
+!               if (lambda_mrcc(i_state,idx_alpha(l_sd)) /= 0.d0) then
+!                 loop = .False.
+!                 exit
+!               endif
+!             enddo
+            loop = .false.
+            if (.not.loop) then
+              call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
+              hIl = hij_mrcc(idx_alpha(l_sd),i_I)
+!                             call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
+              do i_state=1,N_states
+                dka(i_state) = dij(i_I, idx_alpha(l_sd), i_state) * phase * phase2
+                !dka(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(l_sd)), N_int) * phase * phase2 !hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
+                !dka(i_state) = psi_non_ref_coef(idx_alpha(l_sd), i_state) / psi_ref_coef(i_I, i_state) * phase * phase2 
+              enddo
+            endif
+
+            exit
+          endif
+        enddo
+        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)
+!        call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla)
+        do i_state=1,N_states
+          dIa_hla(i_state,k_sd) = dIa(i_state) * hla
+        enddo
+      enddo
+      call omp_set_lock( psi_ref_lock(i_I) )
+      do i_state=1,N_states
+        if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
+          do l_sd=1,idx_alpha(0)
+            k_sd = idx_alpha(l_sd)
+            delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(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)
+          enddo
+        else
+          delta_ii_(i_state,i_I)  = 0.d0
+          do l_sd=1,idx_alpha(0)
+            k_sd = idx_alpha(l_sd)
+            delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + 0.5d0*dIa_hla(i_state,k_sd)
+          enddo
+        endif
+      enddo
+      call omp_unset_lock( psi_ref_lock(i_I) )
+    enddo
+  enddo
+  deallocate (dIa_hla,hij_cache)
+  deallocate(miniList, idx_miniList)
+end
+
+
+
+
+ BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ii, (N_states, N_det_ref) ]
+  use bitmasks
+  implicit none
+  integer :: i, j, i_state
+  
+  !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
+  
+  do i_state = 1, N_states
+    if(mrmode == 3) then
+    do i = 1, N_det_ref
+      delta_ii(i_state,i)= delta_ii_mrcc(i_state,i)
+      do j = 1, N_det_non_ref
+        delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i)
+      end do
+    end do
+!       
+!       do i = 1, N_det_ref
+!         delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state)
+!         do j = 1, N_det_non_ref
+!           delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state)
+!         end do
+!       end do
+    else if(mrmode == 2) then
+      do i = 1, N_det_ref
+        delta_ii(i_state,i)= delta_ii_old(i_state,i)
+        do j = 1, N_det_non_ref
+          delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i)
+        end do
+      end do
+    else if(mrmode == 1) then
+      do i = 1, N_det_ref
+        delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state)
+        do j = 1, N_det_non_ref
+          delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state)
+        end do
+      end do
+    else
+      stop "invalid mrmode"
+    end if
+  end do
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, HP, (2,N_det_non_ref) ]
+  integer :: i
+  do i=1,N_det_non_ref
+    call getHP(psi_non_ref(1,1,i), HP(1,i), HP(2,i), N_int)
+  end do 
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, cepa0_shortcut, (0:N_det_non_ref+1) ]
+&BEGIN_PROVIDER [ integer, det_cepa0_idx, (N_det_non_ref) ]
+&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0_active, (N_int,2,N_det_non_ref) ]
+&BEGIN_PROVIDER [ integer(bit_kind), det_ref_active, (N_int,2,N_det_ref) ]
+&BEGIN_PROVIDER [ integer(bit_kind), active_sorb, (N_int,2) ]
+&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0, (N_int,2,N_det_non_ref) ]
+&BEGIN_PROVIDER [ integer, nlink, (N_det_ref) ]
+&BEGIN_PROVIDER [ integer, linked, (N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ integer, blokMwen, (N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, searchance, (N_det_ref) ]
+&BEGIN_PROVIDER [ integer, child_num, (N_det_non_ref,N_det_ref) ]
+
+  use bitmasks
+  implicit none
+  
+  integer(bit_kind),allocatable :: det_noactive(:,:,:)
+  integer, allocatable :: shortcut(:), idx(:)
+  integer(bit_kind) :: nonactive_sorb(N_int,2), det(N_int, 2)
+  integer i, II, j, k, n, ni, blok, degree
+  logical, external :: detEq
+  
+  allocate(det_noactive(N_int, 2, N_det_non_ref))
+  allocate(idx(N_det_non_ref), shortcut(0:N_det_non_ref+1))
+  print *, "pre start"
+  active_sorb(:,:) = 0_8
+  nonactive_sorb(:,:) = not(0_8)
+  
+  if(N_det_ref > 1) then
+    do i=1, N_det_ref
+    do k=1, N_int
+      active_sorb(k,1) = ior(psi_ref(k,1,i), active_sorb(k,1))
+      active_sorb(k,2) = ior(psi_ref(k,2,i), active_sorb(k,2))
+      nonactive_sorb(k,1) = iand(psi_ref(k,1,i), nonactive_sorb(k,1))
+      nonactive_sorb(k,2) = iand(psi_ref(k,2,i), nonactive_sorb(k,2))
+    end do
+    end do
+    do k=1, N_int
+      active_sorb(k,1) = iand(active_sorb(k,1), not(nonactive_sorb(k,1)))
+      active_sorb(k,2) = iand(active_sorb(k,2), not(nonactive_sorb(k,2)))
+    end do
+  end if
+  
+  
+  do i=1, N_det_non_ref
+  do k=1, N_int
+    det_noactive(k,1,i) = iand(psi_non_ref(k,1,i), not(active_sorb(k,1)))
+    det_noactive(k,2,i) = iand(psi_non_ref(k,2,i), not(active_sorb(k,2))) 
+  end do
+  end do
+  
+  call sort_dets_ab(det_noactive, det_cepa0_idx, cepa0_shortcut, N_det_non_ref, N_int)
+  
+  do i=1,N_det_non_ref
+    det_cepa0(:,:,i) = psi_non_ref(:,:,det_cepa0_idx(i))
+  end do
+  
+  cepa0_shortcut(0) = 1
+  cepa0_shortcut(1) = 1
+  do i=2,N_det_non_ref
+    if(.not. detEq(det_noactive(1,1,i), det_noactive(1,1,i-1), N_int)) then
+      cepa0_shortcut(0) += 1
+      cepa0_shortcut(cepa0_shortcut(0)) = i
+    end if
+  end do
+  cepa0_shortcut(cepa0_shortcut(0)+1) = N_det_non_ref+1
+  
+  if(.true.) then 
+  do i=1,cepa0_shortcut(0)
+    n = cepa0_shortcut(i+1) - cepa0_shortcut(i)
+    call sort_dets_ab(det_cepa0(1,1,cepa0_shortcut(i)), idx, shortcut, n, N_int)
+    do k=1,n
+      idx(k) = det_cepa0_idx(cepa0_shortcut(i)-1+idx(k))
+    end do
+    det_cepa0_idx(cepa0_shortcut(i):cepa0_shortcut(i)+n-1) = idx(:n)
+  end do
+  end if
+  
+  
+  do i=1,N_det_ref
+  do k=1, N_int
+    det_ref_active(k,1,i) = iand(psi_ref(k,1,i), active_sorb(k,1))
+    det_ref_active(k,2,i) = iand(psi_ref(k,2,i), active_sorb(k,2))
+  end do
+  end do
+  
+  do i=1,N_det_non_ref
+  do k=1, N_int
+    det_cepa0_active(k,1,i) = iand(psi_non_ref(k,1,det_cepa0_idx(i)), active_sorb(k,1))
+    det_cepa0_active(k,2,i) = iand(psi_non_ref(k,2,det_cepa0_idx(i)), active_sorb(k,2))
+  end do
+  end do
+
+  do i=1,N_det_non_ref
+    if(.not. detEq(psi_non_ref(1,1,det_cepa0_idx(i)), det_cepa0(1,1,i),N_int)) stop "STOOOP"
+  end do
+  
+  searchance = 0d0
+  child_num = 0
+  do J = 1, N_det_ref
+    nlink(J) = 0
+    do blok=1,cepa0_shortcut(0)
+    do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
+      call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int)
+      if(degree <= 2) then
+        nlink(J) += 1
+        linked(nlink(J),J) = k
+        child_num(k, J) = nlink(J)
+        blokMwen(nlink(J),J) = blok
+        searchance(J) += 1d0 + log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok)))
+      end if
+    end do
+    end do
+  end do
+  print *, "pre done"
+END_PROVIDER
+
+ 
+!  BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ]
+!   use bitmasks
+!   implicit none
+!   integer :: i,j,k
+!   double precision :: Hjk, Hki, Hij, pre(N_det_ref), wall
+!   integer :: i_state, degree, npre, ipre(N_det_ref), npres(N_det_ref)
+!   
+! !   provide lambda_mrcc
+!   npres = 0
+!   delta_cas = 0d0
+!   call wall_time(wall)
+!   print *, "dcas ", wall
+!   do i_state = 1, N_states
+!     !!$OMP PARALLEL DO default(none) schedule(dynamic) private(pre,npre,ipre,j,k,Hjk,Hki,degree) shared(npres,lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref)
+!     do k=1,N_det_non_ref
+!       if(lambda_mrcc(i_state, k) == 0d0) cycle
+!       npre = 0
+!       do i=1,N_det_ref
+!         call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki)
+!         if(Hki /= 0d0) then
+!           !!$OMP ATOMIC
+!           npres(i) += 1
+!           npre += 1
+!           ipre(npre) = i
+!           pre(npre) = Hki
+!         end if
+!       end do
+!       
+!       
+!       do i=1,npre
+!       do j=1,i
+!         !!$OMP ATOMIC
+!         delta_cas(ipre(i),ipre(j),i_state) += pre(i) * pre(j) * lambda_mrcc(i_state, k)
+!       end do
+!       end do
+!     end do
+!     !!$OMP END PARALLEL DO
+!     npre=0
+!     do i=1,N_det_ref
+!       npre += npres(i)
+!     end do
+!     !stop
+!     do i=1,N_det_ref
+!     do j=1,i
+!         delta_cas(j,i,i_state) = delta_cas(i,j,i_state)
+!     end do
+!     end do 
+!   end do
+!   
+!   call wall_time(wall)
+!   print *, "dcas", wall
+! !   stop
+!  END_PROVIDER
+ 
+ 
+ BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ]
+  use bitmasks
+  implicit none
+  integer :: i,j,k
+  double precision :: Hjk, Hki, Hij
+  !double precision, external :: get_dij
+  integer i_state, degree
+   
+  provide lambda_mrcc dIj
+  do i_state = 1, N_states
+    !$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref,dij)
+    do i=1,N_det_ref
+      do j=1,i
+        call get_excitation_degree(psi_ref(1,1,i), psi_ref(1,1,j), degree, N_int)
+        delta_cas(i,j,i_state) = 0d0
+        do k=1,N_det_non_ref
+
+          call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk)
+          
+          delta_cas(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k)
+          !print *, Hjk * get_dij(psi_ref(1,1,i), psi_non_ref(1,1,k), N_int), Hki * get_dij(psi_ref(1,1,j), psi_non_ref(1,1,k), N_int)
+        end do
+        delta_cas(j,i,i_state) = delta_cas(i,j,i_state)
+      end do
+    end do
+    !$OMP END PARALLEL DO
+  end do
+ END_PROVIDER
+ 
+ 
+
+
+logical function isInCassd(a,Nint)
+  use bitmasks
+  implicit none
+  integer, intent(in) :: Nint
+  integer(bit_kind), intent(in) :: a(Nint,2)
+  integer(bit_kind) :: inac, virt
+  integer :: ni, i, deg
+
+  
+  isInCassd = .false.
+     
+  deg = 0
+   do i=1,2
+    do ni=1,Nint
+      virt = iand(not(HF_bitmask(ni,i)), not(active_sorb(ni,i)))
+      deg += popcnt(iand(virt, a(ni,i)))
+      if(deg > 2) return
+    end do
+  end do
+  
+  deg = 0
+  do i=1,2
+    do ni=1,Nint
+      inac = iand(HF_bitmask(ni,i), not(active_sorb(ni,i)))
+      deg += popcnt(xor(iand(inac,a(ni,i)), inac))
+      if(deg > 2) return
+    end do
+  end do
+  isInCassd = .true.
+end function
+
+
+subroutine getHP(a,h,p,Nint)
+  use bitmasks
+  implicit none
+  integer, intent(in) :: Nint
+  integer(bit_kind), intent(in) :: a(Nint,2)
+  integer, intent(out) :: h, p
+  integer(bit_kind) :: inac, virt
+  integer :: ni, i, deg
+
+  
+  !isInCassd = .false.
+  h = 0
+  p = 0
+   
+  deg = 0
+  lp : do i=1,2
+    do ni=1,Nint
+      virt = iand(not(HF_bitmask(ni,i)), not(active_sorb(ni,i)))
+      deg += popcnt(iand(virt, a(ni,i)))
+      if(deg > 2) exit lp
+    end do
+  end do lp
+  p = deg
+
+  deg = 0
+  lh : do i=1,2
+    do ni=1,Nint
+      inac = iand(HF_bitmask(ni,i), not(active_sorb(ni,i)))
+      deg += popcnt(xor(iand(inac,a(ni,i)), inac))
+      if(deg > 2) exit lh
+    end do
+  end do lh
+  h = deg
+  !isInCassd = .true.
+end function
+
+
+ BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ]
+&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii, (N_det_ref,N_states) ]
+  use bitmasks
+  implicit none
+  
+  integer                         :: i_state, i, i_I, J, k, degree, degree2, m, l, deg, ni
+  integer                         :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_, sortRefIdx(N_det_ref)
+  logical                         :: ok
+  double precision                :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1)
+  double precision                :: contrib, contrib2,  HIIi, HJk, wall
+  integer, dimension(0:2,2,2)     :: exc_iI, exc_Ik, exc_IJ
+  integer(bit_kind)               :: det_tmp(N_int, 2), made_hole(N_int,2), made_particle(N_int,2), myActive(N_int,2)
+  integer(bit_kind),allocatable   :: sortRef(:,:,:)
+  integer, allocatable            :: idx_sorted_bit(:)
+  integer, external               :: get_index_in_psi_det_sorted_bit, searchDet
+  logical, external               :: is_in_wavefunction, detEq
+  !double precision, external      :: get_dij
+  integer :: II, blok
+  integer*8, save :: notf = 0
+
+  call wall_time(wall)
+  allocate(idx_sorted_bit(N_det), sortRef(N_int,2,N_det_ref))
+  
+  sortRef(:,:,:) = det_ref_active(:,:,:)
+  call sort_det(sortRef, sortRefIdx, N_det_ref, N_int)
+  
+  idx_sorted_bit(:) = -1
+  do i=1,N_det_non_ref
+    idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
+  enddo
+    
+  ! To provide everything
+  contrib = dij(1, 1, 1)
+  
+  do i_state = 1, N_states
+    delta_mrcepa0_ii(:,:) = 0d0
+    delta_mrcepa0_ij(:,:,:) = 0d0
+
+    !$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii)       &
+    !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib2)       &
+    !$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active)     &
+    !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas) &
+    !$OMP shared(notf,i_state, sortRef, sortRefIdx, dij)
+    do blok=1,cepa0_shortcut(0)
+    do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
+      do II=1,N_det_ref
+        call get_excitation_degree(psi_ref(1,1,II),psi_non_ref(1,1,det_cepa0_idx(i)),degree,N_int)
+        if (degree > 2 ) cycle
+        
+        do ni=1,N_int
+          made_hole(ni,1) = iand(det_ref_active(ni,1,II), xor(det_cepa0_active(ni,1,i), det_ref_active(ni,1,II)))
+          made_hole(ni,2) = iand(det_ref_active(ni,2,II), xor(det_cepa0_active(ni,2,i), det_ref_active(ni,2,II)))
+          
+          made_particle(ni,1) = iand(det_cepa0_active(ni,1,i), xor(det_cepa0_active(ni,1,i), det_ref_active(ni,1,II)))
+          made_particle(ni,2) = iand(det_cepa0_active(ni,2,i), xor(det_cepa0_active(ni,2,i), det_ref_active(ni,2,II)))
+        end do
+      
+        
+  kloop: do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 !i
+          !if(lambda_mrcc(i_state, det_cepa0_idx(k)) == 0d0) cycle
+          
+          do ni=1,N_int
+            if(iand(made_hole(ni,1), det_cepa0_active(ni,1,k)) /= 0) cycle kloop
+            if(iand(made_particle(ni,1), det_cepa0_active(ni,1,k)) /= made_particle(ni,1)) cycle kloop
+            if(iand(made_hole(ni,2), det_cepa0_active(ni,2,k)) /= 0) cycle kloop
+            if(iand(made_particle(ni,2), det_cepa0_active(ni,2,k)) /= made_particle(ni,2)) cycle kloop
+          end do
+          do ni=1,N_int
+            myActive(ni,1) = xor(det_cepa0_active(ni,1,k), made_hole(ni,1))
+            myActive(ni,1) = xor(myActive(ni,1), made_particle(ni,1))
+            myActive(ni,2) = xor(det_cepa0_active(ni,2,k), made_hole(ni,2))
+            myActive(ni,2) = xor(myActive(ni,2), made_particle(ni,2))
+          end do
+          
+          j = searchDet(sortRef, myActive, N_det_ref, N_int)
+          if(j == -1) then
+           cycle
+          end if
+          j = sortRefIdx(j)
+          !$OMP ATOMIC
+          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)
+          
+          if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then
+            contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state)
+            !$OMP ATOMIC
+            delta_mrcepa0_ii(J,i_state) -= contrib2 
+          else
+            contrib = contrib * 0.5d0
+          end if
+          !$OMP ATOMIC
+          delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib
+
+        end do kloop
+      end do
+    end do
+    end do
+    !$OMP END PARALLEL DO
+  end do
+  deallocate(idx_sorted_bit)
+  call wall_time(wall)
+  print *, "cepa0", wall, notf
+  !stop
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, delta_sub_ij, (N_det_ref,N_det_non_ref,N_states) ]
+&BEGIN_PROVIDER [ double precision, delta_sub_ii, (N_det_ref, N_states) ]
+  use bitmasks
+  implicit none
+  
+  integer                         :: i_state, i, i_I, J, k, degree, degree2, l, deg, ni
+  integer                         :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_
+  logical                         :: ok
+  double precision                :: phase_Ji, phase_Ik, phase_Ii
+  double precision                :: contrib, contrib2, delta_IJk, HJk, HIk, HIl
+  integer, dimension(0:2,2,2)     :: exc_Ik, exc_Ji, exc_Ii
+  integer(bit_kind)               :: det_tmp(N_int, 2), det_tmp2(N_int, 2)
+  integer, allocatable            :: idx_sorted_bit(:)
+  integer, external               :: get_index_in_psi_det_sorted_bit
+  
+  integer :: II, blok
+  
+  provide delta_cas lambda_mrcc
+  allocate(idx_sorted_bit(N_det))
+  idx_sorted_bit(:) = -1
+  do i=1,N_det_non_ref
+    idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
+  enddo
+    
+  do i_state = 1, N_states
+    delta_sub_ij(:,:,:) = 0d0
+    delta_sub_ii(:,:) = 0d0
+    
+    provide mo_bielec_integrals_in_map
+    
+    
+    !$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_sub_ij, delta_sub_ii)       &
+    !$OMP private(i, J, k, degree, degree2, l, deg, ni)       &
+    !$OMP private(p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_)     &
+    !$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib2, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) &
+    !$OMP private(det_tmp, det_tmp2, II, blok)    &
+    !$OMP shared(idx_sorted_bit, N_det_non_ref, N_det_ref, N_int, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef)   &
+    !$OMP shared(i_state,lambda_mrcc, hf_bitmask, active_sorb)
+    do i=1,N_det_non_ref
+      if(mod(i,1000) == 0) print *, i, "/", N_det_non_ref
+      do J=1,N_det_ref
+        call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_Ji,degree,phase_Ji,N_int)
+        if(degree == -1) cycle
+        
+        
+        do II=1,N_det_ref
+          call apply_excitation(psi_ref(1,1,II),exc_Ji,det_tmp,ok,N_int)
+
+          if(.not. ok) cycle
+          l = get_index_in_psi_det_sorted_bit(det_tmp, N_int)
+          if(l == 0) cycle
+          l = idx_sorted_bit(l)
+          
+          call i_h_j(psi_ref(1,1,II), det_tmp, N_int, HIl)
+          
+          do k=1,N_det_non_ref
+            if(lambda_mrcc(i_state, k) == 0d0) cycle
+            call get_excitation(psi_ref(1,1,II),psi_non_ref(1,1,k),exc_Ik,degree2,phase_Ik,N_int)
+            
+            det_tmp(:,:) = 0_bit_kind
+            det_tmp2(:,:) = 0_bit_kind
+            
+            ok = .true.
+            do ni=1,N_int
+              det_tmp(ni,1) = iand(xor(HF_bitmask(ni,1), psi_non_ref(ni,1,k)), not(active_sorb(ni,1)))
+              det_tmp(ni,2) = iand(xor(HF_bitmask(ni,1), psi_non_ref(ni,1,i)), not(active_sorb(ni,1)))
+              ok = ok .and. (popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) == popcnt(xor(det_tmp(ni,1), det_tmp(ni,2))))
+
+              det_tmp(ni,1) = iand(xor(HF_bitmask(ni,2), psi_non_ref(ni,2,k)), not(active_sorb(ni,2)))
+              det_tmp(ni,2) = iand(xor(HF_bitmask(ni,2), psi_non_ref(ni,2,i)), not(active_sorb(ni,2)))
+              ok = ok .and. (popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) == popcnt(xor(det_tmp(ni,1), det_tmp(ni,2))))
+            end do
+            
+            if(ok) cycle
+            
+            
+            call i_h_j(psi_ref(1,1,J), psi_non_ref(1,1,k), N_int, HJk)
+            call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,k), N_int, HIk)
+            if(HJk == 0) cycle
+            !assert HIk == 0
+            delta_IJk = HJk * HIk * lambda_mrcc(i_state, k)
+            call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
+            if(ok) cycle
+            contrib = delta_IJk * HIl * lambda_mrcc(i_state,l)   
+            if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then
+              contrib2 = contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state)
+              !$OMP ATOMIC
+              delta_sub_ii(II,i_state) -= contrib2
+            else
+              contrib = contrib * 0.5d0
+            endif
+            !$OMP ATOMIC
+            delta_sub_ij(II, i, i_state) += contrib
+          end do
+        end do
+      end do
+    end do
+    !$OMP END PARALLEL DO
+  end do
+  deallocate(idx_sorted_bit)
+END_PROVIDER
+
+
+subroutine set_det_bit(det, p, s)
+   implicit none
+   integer(bit_kind),intent(inout) :: det(N_int, 2)
+   integer, intent(in) :: p, s
+   integer :: ni, pos
+ 
+   ni = (p-1)/bit_kind_size + 1
+   pos = mod(p-1, bit_kind_size)
+   det(ni,s) = ibset(det(ni,s), pos)
+end subroutine
+
+
+BEGIN_PROVIDER [ double precision, h_, (N_det_ref,N_det_non_ref) ]
+  implicit none
+  integer :: i,j
+  do i=1,N_det_ref
+  do j=1,N_det_non_ref
+    call i_h_j(psi_ref(1,1,i), psi_non_ref(1,1,j), N_int, h_(i,j))
+  end do
+  end do
+END_PROVIDER
+
+
+
+subroutine filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList)
+
+ use bitmasks
+ implicit none
+
+  integer, intent(in)            :: i_generator,n_selected, Nint
+
+  integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
+  integer                        :: i,j,k,m
+  logical                        :: is_in_wavefunction
+  integer,allocatable            :: degree(:)
+  integer,allocatable            :: idx(:)
+  logical                        :: good
+
+  integer(bit_kind), intent(inout) :: tq(Nint,2,n_selected) !! intent(out)
+  integer, intent(out)           :: N_tq
+  
+  integer                        :: nt,ni
+  logical, external              :: is_connected_to, is_generable
+  
+  integer(bit_kind),intent(in)  :: miniList(Nint,2,N_det_generators)
+  integer,intent(in)            :: N_miniList
+  
+  allocate(degree(psi_det_size)) 
+  allocate(idx(0:psi_det_size))
+  N_tq = 0
+
+  i_loop : do i=1,N_selected
+    do k=1, N_minilist
+      if(is_generable(miniList(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop
+    end do
+
+    ! Select determinants that are triple or quadruple excitations
+    ! from the ref
+    good = .True.
+    call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) 
+    !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector
+    do k=1,idx(0)
+      if (degree(k) < 3) then
+        good = .False.
+        exit
+      endif
+    enddo
+    if (good) then
+      if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then
+        N_tq += 1
+        do k=1,N_int
+          tq(k,1,N_tq) = det_buffer(k,1,i)
+          tq(k,2,N_tq) = det_buffer(k,2,i)
+        enddo
+      endif
+    endif
+  enddo i_loop
+end
+
+
+subroutine filter_tq_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
+
+ use bitmasks
+ implicit none
+
+  integer, intent(in)            :: i_generator,n_selected, Nint
+
+  integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
+  integer                        :: i,j,k,m
+  logical                        :: is_in_wavefunction
+  integer,allocatable            :: degree(:)
+  integer,allocatable            :: idx(:)
+  logical                        :: good
+
+  integer(bit_kind), intent(inout) :: tq(Nint,2,n_selected) !! intent(out)
+  integer, intent(out)           :: N_tq
+  
+  integer                        :: nt,ni
+  logical, external              :: is_connected_to, is_generable
+
+  integer(bit_kind),intent(in)  :: microlist(Nint,2,*)
+  integer,intent(in)  :: ptr_microlist(0:*)
+  integer,intent(in)            :: N_microlist(0:*)
+  integer(bit_kind),intent(in)   :: key_mask(Nint, 2)
+  
+  integer :: mobiles(2), smallerlist
+  
+  
+  allocate(degree(psi_det_size)) 
+  allocate(idx(0:psi_det_size))
+  N_tq = 0
+  
+  i_loop : do i=1,N_selected
+    call getMobiles(det_buffer(1,1,i), key_mask, mobiles, Nint) 
+    if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then
+      smallerlist = mobiles(1)
+    else
+      smallerlist = mobiles(2)
+    end if
+    
+    if(N_microlist(smallerlist) > 0) then
+      do k=ptr_microlist(smallerlist), ptr_microlist(smallerlist)+N_microlist(smallerlist)-1
+        if(is_generable(microlist(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop
+      end do
+    end if
+    
+    if(N_microlist(0) > 0) then
+      do k=1, N_microlist(0)
+        if(is_generable(microlist(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop
+      end do
+    end if
+
+    ! Select determinants that are triple or quadruple excitations
+    ! from the ref
+    good = .True.
+    call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) 
+    !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector
+    do k=1,idx(0)
+      if (degree(k) < 3) then
+        good = .False.
+        exit
+      endif
+    enddo
+    if (good) then
+      if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then
+        N_tq += 1
+        do k=1,N_int
+          tq(k,1,N_tq) = det_buffer(k,1,i)
+          tq(k,2,N_tq) = det_buffer(k,2,i)
+        enddo
+      endif
+    endif
+  enddo i_loop
+end
+
+
+
+
diff --git a/plugins/mrcc_selected/dressing_slave.irp.f b/plugins/mrcc_selected/dressing_slave.irp.f
new file mode 100644
index 00000000..f1d6f029
--- /dev/null
+++ b/plugins/mrcc_selected/dressing_slave.irp.f
@@ -0,0 +1,593 @@
+subroutine mrsc2_dressing_slave_tcp(i)
+  implicit none
+  integer, intent(in)            :: i
+  BEGIN_DOC
+! Task for parallel MR-SC2
+  END_DOC
+  call mrsc2_dressing_slave(0,i)
+end
+
+
+subroutine mrsc2_dressing_slave_inproc(i)
+  implicit none
+  integer, intent(in)            :: i
+  BEGIN_DOC
+! Task for parallel MR-SC2
+  END_DOC
+  call mrsc2_dressing_slave(1,i)
+end
+
+subroutine mrsc2_dressing_slave(thread,iproc)
+  use f77_zmq
+
+  implicit none
+  BEGIN_DOC
+! Task for parallel MR-SC2
+  END_DOC
+  integer,  intent(in)            :: thread, iproc
+!   integer                        :: j,l
+  integer                        :: rc 
+
+  integer                        :: worker_id, task_id 
+  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 
+
+  double precision, allocatable  :: delta(:,:,:) 
+  
+
+
+  integer                         :: i_state, i, i_I, J, k, k2, k1, kk, ll, degree, degree2, 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_iI, phase_Ik, phase_Jl, phase_Ji, phase_al
+  double precision                :: diI, hIi, hJi, delta_JI, dkI, HkI, ci_inv(N_states), cj_inv(N_states)
+  double precision                :: contrib, wall, iwall
+  double precision, allocatable   :: dleat(:,:,:)
+  integer, dimension(0:2,2,2)     :: exc_iI, exc_Ik, exc_IJ
+  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, isInCassd, detEq
+  integer,allocatable :: komon(:)
+  logical :: komoned
+  !double precision, external :: get_dij
+     
+  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)
+
+  allocate (dleat(N_states, N_det_non_ref, 2), delta(N_states,0:N_det_non_ref, 2))
+  allocate(komon(0:N_det_non_ref))
+
+  do 
+    call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task)
+    if (task_id == 0) exit
+    read (task,*) i_I, J, k1, k2
+    do i_state=1, N_states
+      ci_inv(i_state) = 1.d0 / psi_ref_coef(i_I,i_state)
+      cj_inv(i_state) = 1.d0 / psi_ref_coef(J,i_state)
+    end do
+    !delta = 0.d0
+    n = 0
+    delta(:,0,:) = 0d0
+    delta(:,:nlink(J),1) = 0d0
+    delta(:,:nlink(i_I),2) = 0d0
+    komon(0) = 0
+    komoned = .false.
+    
+    
+    
+    
+    do kk = k1, k2
+      k = det_cepa0_idx(linked(kk, i_I))
+      blok = blokMwen(kk, i_I)
+      
+      call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int)
+    
+      if(J /= i_I) then
+        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)
+      else
+        l = k
+        ll = kk
+      end if
+      
+        
+      if(.not. komoned) then
+        m = 0
+        m2 = 0
+          
+        do while(m < nlink(i_I) .and. m2 < nlink(J))
+          m += 1
+          m2 += 1
+          if(linked(m, i_I) < linked(m2, J)) then
+            m2 -= 1
+            cycle
+          else if(linked(m, i_I) > linked(m2, J)) then
+            m -= 1
+            cycle
+          end if
+          i = det_cepa0_idx(linked(m, i_I))
+          
+          if(h_(J,i) == 0.d0) cycle
+          if(h_(i_I,i) == 0.d0) cycle
+          
+          !ok = .false.
+          !do i_state=1, N_states
+          !  if(lambda_mrcc(i_state, i) /= 0d0) then
+          !    ok = .true.
+          !    exit
+          !  end if
+          !end do
+          !if(.not. ok) cycle
+!         
+          
+          komon(0) += 1
+          kn = komon(0)
+          komon(kn) = i
+          
+          
+!           call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ji,N_int)
+!           if(I_i /= J) call get_excitation(psi_ref(1,1,I_i),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ii,N_int)
+!           if(I_i == J) phase_Ii = phase_Ji
+          
+          do i_state = 1,N_states
+            dkI = h_(J,i) * dij(i_I, i, i_state)!get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,i), N_int)
+            !dkI = h_(J,i) * h_(i_I,i) * lambda_mrcc(i_state, i)
+            dleat(i_state, kn, 1) = dkI
+            dleat(i_state, kn, 2) = dkI
+          end do
+
+        end do
+          
+        komoned = .true.
+      end if
+      
+      
+      do m = 1, komon(0)
+        
+        i = komon(m)
+        
+        call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
+        if(.not. ok) cycle
+        if(HP(1,i) + HP(1,k) <= 2 .and. HP(2,i) + HP(2,k) <= 2) then
+!           if(is_in_wavefunction(det_tmp, N_int)) cycle
+          cycle
+        end if
+        
+        !if(isInCassd(det_tmp, N_int)) cycle
+          
+        do i_state = 1, N_states 
+          !if(lambda_mrcc(i_state, i) == 0d0) cycle
+          
+
+          !contrib = h_(i_I,k) * lambda_mrcc(i_state, k) * dleat(i_state, m, 2)! * phase_al
+          contrib =  dij(i_I, k, i_state) * dleat(i_state, m, 2)
+          delta(i_state,ll,1) += contrib
+          if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then
+            delta(i_state,0,1) -= contrib * ci_inv(i_state) * psi_non_ref_coef(l,i_state)
+          endif
+          
+          if(I_i == J) cycle
+          !contrib = h_(J,l) * lambda_mrcc(i_state, l) * dleat(i_state, m, 1)! * phase_al
+          contrib =  dij(J, l, i_state) * dleat(i_state, m, 1)
+          delta(i_state,kk,2) += contrib
+          if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then
+            delta(i_state,0,2) -= contrib * cj_inv(i_state) * psi_non_ref_coef(k,i_state)
+          end if
+        enddo !i_state
+      end do ! while
+    end do ! kk
+
+      
+      call push_mrsc2_results(zmq_socket_push, I_i, J, delta, task_id) 
+      call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
+ 
+!     end if
+    
+  enddo
+
+  deallocate(delta)
+
+  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)
+
+end
+
+
+subroutine push_mrsc2_results(zmq_socket_push, I_i, J, delta, task_id) 
+  use f77_zmq 
+  implicit none 
+  BEGIN_DOC 
+! Push integrals in the push socket 
+  END_DOC 
+
+  integer, intent(in)            :: i_I, J
+  integer(ZMQ_PTR), intent(in)   :: zmq_socket_push 
+  double precision,intent(inout) :: delta(N_states, 0:N_det_non_ref, 2)
+  integer, intent(in)            :: task_id 
+  integer                        :: rc , i_state,  i, kk, li
+  integer,allocatable            :: idx(:,:)
+  integer                        :: n(2)
+  logical :: ok
+  
+  allocate(idx(N_det_non_ref,2))
+  rc = f77_zmq_send( zmq_socket_push, i_I, 4, ZMQ_SNDMORE) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_send( zmq_socket_push, i_I, 4, ZMQ_SNDMORE)' 
+    stop 'error' 
+  endif 
+
+  rc = f77_zmq_send( zmq_socket_push, J, 4, ZMQ_SNDMORE) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_send( zmq_socket_push, J, 4, ZMQ_SNDMORE)' 
+    stop 'error' 
+  endif 
+  
+  
+  do kk=1,2
+    n(kk)=0
+    if(kk == 1) li = nlink(j)
+    if(kk == 2) li = nlink(i_I)
+    do i=1, li
+      ok = .false.
+      do i_state=1,N_states
+        if(delta(i_state, i, kk) /= 0d0) then
+          ok = .true.
+          exit
+        end if
+      end do
+      
+      if(ok) then
+        n(kk) += 1
+!         idx(n,kk) = i
+        if(kk == 1) then
+          idx(n(1),1) = det_cepa0_idx(linked(i, J))
+        else
+          idx(n(2),2) = det_cepa0_idx(linked(i, i_I))
+        end if
+        
+        do i_state=1, N_states
+            delta(i_state, n(kk), kk) = delta(i_state, i, kk)
+        end do
+      end if
+    end do
+    
+    rc = f77_zmq_send( zmq_socket_push, n(kk), 4, ZMQ_SNDMORE) 
+    if (rc /= 4) then 
+      print *, irp_here,  'f77_zmq_send( zmq_socket_push, n, 4, ZMQ_SNDMORE)' 
+      stop 'error' 
+    endif 
+    
+    if(n(kk) /= 0) then
+      rc = f77_zmq_send( zmq_socket_push, delta(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE)  ! delta(1,0,1) = delta_I   delta(1,0,2) = delta_J 
+      if (rc /=  (n(kk)+1)*8*N_states) then 
+        print *, irp_here,  'f77_zmq_send( zmq_socket_push, delta, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)'
+        stop 'error' 
+      endif 
+      
+      rc = f77_zmq_send( zmq_socket_push, idx(1,kk), n(kk)*4, ZMQ_SNDMORE) 
+      if (rc /=  n(kk)*4) then 
+        print *, irp_here,  'f77_zmq_send( zmq_socket_push, delta, 8*n(kk), ZMQ_SNDMORE)'
+        stop 'error' 
+      endif 
+    end if
+  end do
+  
+  
+  rc = f77_zmq_send( zmq_socket_push, task_id, 4, 0) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_send( zmq_socket_push, task_id, 4, 0)' 
+    stop 'error' 
+  endif 
+ 
+! ! Activate is zmq_socket_push is a REQ 
+!   integer :: idummy 
+!   rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0) 
+!   if (rc /= 4) then 
+!     print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)' 
+!     stop 'error' 
+!   endif 
+end
+
+
+
+subroutine pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, task_id) 
+  use f77_zmq 
+  implicit none 
+  BEGIN_DOC 
+! Push integrals in the push socket 
+  END_DOC 
+
+  integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull 
+  integer, intent(out)           :: i_I, J, n(2)
+  double precision, intent(inout) :: delta(N_states, 0:N_det_non_ref, 2)
+  integer, intent(out)           :: task_id 
+  integer                        :: rc , i, kk
+  integer,intent(inout) :: idx(N_det_non_ref,2)
+  logical :: ok
+  
+  rc = f77_zmq_recv( zmq_socket_pull, i_I, 4, ZMQ_SNDMORE) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, i_I, 4, ZMQ_SNDMORE)' 
+    stop 'error' 
+  endif 
+  
+  rc = f77_zmq_recv( zmq_socket_pull, J, 4, ZMQ_SNDMORE) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, J, 4, ZMQ_SNDMORE)' 
+    stop 'error' 
+  endif 
+  
+  do kk = 1, 2
+    rc = f77_zmq_recv( zmq_socket_pull, n(kk), 4, ZMQ_SNDMORE) 
+    if (rc /= 4) then 
+      print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, n, 4, ZMQ_SNDMORE)' 
+      stop 'error' 
+    endif 
+    
+    if(n(kk) /= 0) then
+      rc = f77_zmq_recv( zmq_socket_pull, delta(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE) 
+      if (rc /= (n(kk)+1)*8*N_states) then 
+        print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, delta, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)'
+        stop 'error' 
+      endif 
+      
+      rc = f77_zmq_recv( zmq_socket_pull, idx(1,kk), n(kk)*4, ZMQ_SNDMORE) 
+      if (rc /= n(kk)*4) then 
+        print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, delta, n(kk)*4, ZMQ_SNDMORE)'
+        stop 'error' 
+      endif 
+    end if
+  end do
+  
+  rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)' 
+    stop 'error' 
+  endif 
+
+ 
+! ! Activate is zmq_socket_pull is a REP 
+!   integer :: idummy 
+!   rc = f77_zmq_send( zmq_socket_pull, idummy, 4, 0) 
+!   if (rc /= 4) then 
+!     print *, irp_here, 'f77_zmq_send( zmq_socket_pull, idummy, 4, 0)' 
+!     stop 'error' 
+!   endif 
+end
+
+
+
+subroutine mrsc2_dressing_collector(delta_ii_,delta_ij_)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC 
+! Collects results from the AO integral calculation 
+  END_DOC 
+ 
+  double precision,intent(inout)               :: delta_ij_(N_states,N_det_non_ref,N_det_ref) 
+  double precision,intent(inout)               :: delta_ii_(N_states,N_det_ref)
+
+!   integer                        :: j,l
+  integer                        :: rc 
+   
+  double precision, allocatable  :: delta(:,:,:) 
+   
+  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*8                      :: control, accu 
+  integer                        :: task_id, more 
+ 
+  integer                        :: I_i, J, l, i_state, n(2), kk
+  integer,allocatable :: idx(:,:)
+  
+  delta_ii_(:,:) = 0d0
+  delta_ij_(:,:,:) = 0d0
+
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() 
+  zmq_socket_pull = new_zmq_pull_socket() 
+ 
+  allocate ( delta(N_states,0:N_det_non_ref,2) ) 
+  
+  allocate(idx(N_det_non_ref,2))
+  more = 1 
+  do while (more == 1) 
+         
+    call pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, task_id)
+    
+
+      do l=1, n(1)
+        do i_state=1,N_states
+          delta_ij_(i_state,idx(l,1),i_I) += delta(i_state,l,1)
+        end do
+      end do
+      
+      do l=1, n(2)
+        do i_state=1,N_states
+          delta_ij_(i_state,idx(l,2),J) += delta(i_state,l,2)
+        end do
+      end do
+
+    
+!     
+!       do l=1,nlink(J)
+!         do i_state=1,N_states
+!           delta_ij_(i_state,det_cepa0_idx(linked(l,J)),i_I) += delta(i_state,l,1)
+!           delta_ij_(i_state,det_cepa0_idx(linked(l,i_I)),j) += delta(i_state,l,2)
+!         end do
+!       end do
+!      
+       if(n(1) /= 0) then 
+       do i_state=1,N_states
+         delta_ii_(i_state,i_I) += delta(i_state,0,1)
+       end do
+       end if
+
+      if(n(2) /= 0) then 
+       do i_state=1,N_states
+         delta_ii_(i_state,J) += delta(i_state,0,2)
+       end do
+       end if
+
+
+    if (task_id /= 0) then 
+      call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) 
+    endif 
+
+ 
+  enddo 
+  deallocate( delta )
+ 
+  call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
+  call end_zmq_pull_socket(zmq_socket_pull)
+
+end
+
+
+
+
+ BEGIN_PROVIDER [ double precision, delta_ij_old, (N_states,N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ii_old, (N_states,N_det_ref) ]
+  implicit none
+  
+  integer                         :: i_state, i, i_I, J, k, kk, degree, degree2, m, l, deg, ni, m2
+  integer                         :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s, nex, nzer, ntot
+!   integer, allocatable            :: linked(:,:), blokMwen(:, :), nlink(:)
+  logical                         :: ok
+  double precision                :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al, diI, hIi, hJi, delta_JI, dkI(N_states), HkI, ci_inv(N_states), dia_hla(N_states)
+  double precision                :: contrib, wall, iwall ! , searchance(N_det_ref)
+  integer, dimension(0:2,2,2)     :: exc_iI, exc_Ik, exc_IJ
+  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, isInCassd, detEq
+  character*(512)                :: task 
+  integer(ZMQ_PTR)               :: zmq_to_qp_run_socket 
+  
+  integer :: KKsize = 1000000
+  
+  
+  call new_parallel_job(zmq_to_qp_run_socket,'mrsc2')
+
+
+  call wall_time(iwall)
+!   allocate(linked(N_det_non_ref, N_det_ref), blokMwen(N_det_non_ref, N_det_ref), nlink(N_det_ref))
+  
+
+!   searchance = 0d0
+!   do J = 1, N_det_ref
+!     nlink(J) = 0
+!     do blok=1,cepa0_shortcut(0)
+!     do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
+!       call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int)
+!       if(degree <= 2) then
+!         nlink(J) += 1
+!         linked(nlink(J),J) = k
+!         blokMwen(nlink(J),J) = blok
+!         searchance(J) += 1d0 + log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok)))
+!       end if
+!     end do
+!     end do
+!   end do
+
+     
+     
+!   stop
+  nzer = 0
+  ntot = 0
+  do nex = 3, 0, -1
+    print *, "los ",nex
+    do I_s = N_det_ref, 1, -1
+!         if(mod(I_s,1) == 0) then
+!           call wall_time(wall)
+!           wall = wall-iwall
+!           print *, I_s, "/", N_det_ref, wall * (dfloat(N_det_ref) / dfloat(I_s)), wall, wall * (dfloat(N_det_ref) / dfloat(I_s))-wall
+!         end if
+
+
+      do J_s = 1, I_s
+      
+        call get_excitation_degree(psi_ref(1,1,J_s), psi_ref(1,1,I_s), degree, N_int)
+        if(degree /= nex) cycle
+        if(nex == 3) nzer = nzer + 1
+        ntot += 1
+!           if(degree > 3) then
+!             deg += 1
+!             cycle
+!           else if(degree == -10) then
+!             KKsize = 100000
+!           else
+!             KKsize = 1000000
+!           end if
+        
+        
+        
+        if(searchance(I_s) < searchance(J_s)) then
+          i_I = I_s
+          J = J_s
+        else
+          i_I = J_s
+          J = I_s
+        end if
+        
+        KKsize = nlink(1)
+        if(nex == 0) KKsize = int(float(nlink(1)) / float(nlink(i_I)) * (float(nlink(1)) / 64d0))
+        
+        !if(KKsize == 0) stop "ZZEO"
+        
+        do kk = 1 , nlink(i_I), KKsize
+          write(task,*) I_i, J, kk, int(min(kk+KKsize-1, nlink(i_I)))
+          call add_task_to_taskserver(zmq_to_qp_run_socket,task)
+        end do
+        
+  !       do kk = 1 , nlink(i_I)
+  !         k = linked(kk,i_I)
+  !         blok = blokMwen(kk,i_I)
+  !         write(task,*) I_i, J, k, blok
+  !         call add_task_to_taskserver(zmq_to_qp_run_socket,task)
+  !         
+  !       enddo !kk
+      enddo !J
+          
+    enddo !I
+  end do ! nex
+  print *, "tasked"
+!  integer(ZMQ_PTR)               ∷ collector_thread
+!  external                       ∷ ao_bielec_integrals_in_map_collector
+!  rc = pthread_create(collector_thread, mrsc2_dressing_collector)
+  print *, nzer, ntot, float(nzer) / float(ntot)
+  provide nproc
+  !$OMP PARALLEL DEFAULT(none)  SHARED(delta_ii_old,delta_ij_old)  PRIVATE(i) NUM_THREADS(nproc+1)
+      i = omp_get_thread_num()
+      if (i==0) then
+        call mrsc2_dressing_collector(delta_ii_old,delta_ij_old)
+      else
+        call mrsc2_dressing_slave_inproc(i)
+      endif
+  !$OMP END PARALLEL
+
+!  rc = pthread_join(collector_thread)
+  call end_parallel_job(zmq_to_qp_run_socket, 'mrsc2')
+    
+
+END_PROVIDER
+
+
+
diff --git a/plugins/mrcc_selected/ezfio_interface.irp.f b/plugins/mrcc_selected/ezfio_interface.irp.f
new file mode 100644
index 00000000..062af449
--- /dev/null
+++ b/plugins/mrcc_selected/ezfio_interface.irp.f
@@ -0,0 +1,61 @@
+! DO NOT MODIFY BY HAND
+! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
+! from file /home/scemama/quantum_package/src/mrcc_selected/EZFIO.cfg
+
+
+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_mrcc_selected_thresh_dressed_ci(has)
+  if (has) then
+    call ezfio_get_mrcc_selected_thresh_dressed_ci(thresh_dressed_ci)
+  else
+    print *, 'mrcc_selected/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_mrcc_selected_n_it_max_dressed_ci(has)
+  if (has) then
+    call ezfio_get_mrcc_selected_n_it_max_dressed_ci(n_it_max_dressed_ci)
+  else
+    print *, 'mrcc_selected/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_mrcc_selected_lambda_type(has)
+  if (has) then
+    call ezfio_get_mrcc_selected_lambda_type(lambda_type)
+  else
+    print *, 'mrcc_selected/lambda_type not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER
diff --git a/plugins/mrcc_selected/mrcc_selected.irp.f b/plugins/mrcc_selected/mrcc_selected.irp.f
new file mode 100644
index 00000000..91592e62
--- /dev/null
+++ b/plugins/mrcc_selected/mrcc_selected.irp.f
@@ -0,0 +1,19 @@
+program mrsc2sub
+  implicit none
+  double precision, allocatable  :: energy(:)
+  allocate (energy(N_states))
+  
+  !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
+  mrmode = 3
+  
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call print_cas_coefs
+  call set_generators_bitmasks_as_holes_and_particles
+  call run(N_states,energy)
+  if(do_pt2_end)then
+    call run_pt2(N_states,energy)
+  endif
+  deallocate(energy)
+end
+
diff --git a/plugins/mrcc_selected/mrcepa0_general.irp.f b/plugins/mrcc_selected/mrcepa0_general.irp.f
new file mode 100644
index 00000000..e3a2d1f5
--- /dev/null
+++ b/plugins/mrcc_selected/mrcepa0_general.irp.f
@@ -0,0 +1,245 @@
+
+
+subroutine run(N_st,energy)
+  implicit none
+  
+  integer, intent(in) :: N_st 
+  double precision, intent(out) :: energy(N_st) 
+
+  integer :: i,j
+
+  double precision :: E_new, E_old, delta_e
+  integer :: iteration
+  double precision :: E_past(4)
+  
+  integer :: n_it_mrcc_max
+  double precision :: thresh_mrcc
+  double precision, allocatable :: lambda(:)
+  allocate (lambda(N_states))
+
+  
+  thresh_mrcc = thresh_dressed_ci
+  n_it_mrcc_max = n_it_max_dressed_ci
+
+  if(n_it_mrcc_max == 1) then
+    do j=1,N_states_diag
+      do i=1,N_det
+        psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
+      enddo
+    enddo
+    SOFT_TOUCH psi_coef ci_energy_dressed
+    call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
+    call ezfio_set_mrcepa0_energy(ci_energy_dressed(1))
+    call save_wavefunction
+    energy(:) = ci_energy_dressed(:)
+  else
+    E_new = 0.d0
+    delta_E = 1.d0
+    iteration = 0
+    lambda = 1.d0
+    do while (delta_E > thresh_mrcc)
+      iteration += 1
+      print *,  '===========================' 
+      print *,  'MRCEPA0 Iteration', iteration
+      print *,  '===========================' 
+      print *,  ''
+      E_old = sum(ci_energy_dressed)
+      call write_double(6,ci_energy_dressed(1),"MRCEPA0 energy")
+      call diagonalize_ci_dressed(lambda)
+      E_new = sum(ci_energy_dressed)
+      delta_E = dabs(E_new - E_old)
+      call save_wavefunction
+      call ezfio_set_mrcepa0_energy(ci_energy_dressed(1))
+      if (iteration >= n_it_mrcc_max) then
+        exit
+      endif
+    enddo
+    call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy")
+    energy(:) = ci_energy_dressed(:)
+  endif
+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_old(N_st,energy) 
+  implicit none 
+  integer :: i,j,k 
+  integer, intent(in)          :: N_st 
+  double precision, intent(in) :: energy(N_st) 
+  double precision :: pt2_redundant(N_st), pt2(N_st)
+  double precision :: norm_pert(N_st),H_pert_diag(N_st)
+  
+  pt2_redundant = 0.d0 
+  pt2 = 0d0
+  !if(lambda_mrcc_pt2(0)  == 0) return
+  
+  print*,'Last iteration only to compute the PT2' 
+  
+  print * ,'Computing the redundant PT2 contribution'
+
+  if (mrmode == 1) then
+
+    N_det_generators = lambda_mrcc_kept(0) 
+    N_det_selectors = lambda_mrcc_kept(0) 
+
+    do i=1,N_det_generators
+      j = lambda_mrcc_kept(i)
+      do k=1,N_int
+        psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
+        psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
+        psi_selectors(k,1,i) = psi_non_ref(k,1,j)
+        psi_selectors(k,2,i) = psi_non_ref(k,2,j)
+      enddo
+      do k=1,N_st
+        psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
+        psi_selectors_coef(i,k) = psi_non_ref_coef(j,k)
+      enddo
+    enddo
+
+  else
+
+    N_det_generators = N_det_non_ref
+    N_det_selectors = N_det_non_ref
+
+    do i=1,N_det_generators
+      j = i
+      do k=1,N_int
+        psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
+        psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
+        psi_selectors(k,1,i) = psi_non_ref(k,1,j)
+        psi_selectors(k,2,i) = psi_non_ref(k,2,j)
+      enddo
+      do k=1,N_st
+        psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
+        psi_selectors_coef(i,k) = psi_non_ref_coef(j,k)
+      enddo
+    enddo
+
+  endif
+  
+  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_redundant, norm_pert, H_pert_diag,  N_st) 
+  
+  print * ,'Computing the remaining contribution'
+
+  threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
+  threshold_generators = max(threshold_generators,threshold_generators_pt2)
+
+  N_det_generators = N_det_non_ref + N_det_ref
+  N_det_selectors = N_det_non_ref + N_det_ref
+  
+  psi_det_generators(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref)
+  psi_selectors(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref)
+  psi_coef_generators(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:)
+  psi_selectors_coef(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:)
+  
+  do i=N_det_ref+1,N_det_generators
+    j = i-N_det_ref
+    do k=1,N_int
+      psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
+      psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
+      psi_selectors(k,1,i) = psi_non_ref(k,1,j)
+      psi_selectors(k,2,i) = psi_non_ref(k,2,j)
+    enddo
+    do k=1,N_st
+      psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
+      psi_selectors_coef(i,k) = psi_non_ref_coef(j,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) 
+ 
+
+  print *, "Redundant PT2 :",pt2_redundant
+  print *, "Full      PT2 :",pt2
+  print *, lambda_mrcc_kept(0), N_det, N_det_ref, psi_coef(1,1), psi_ref_coef(1,1)
+  pt2 = pt2 - pt2_redundant
+  
+  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_mrcepa0_energy_pt2(energy(1)+pt2(1))
+
+end 
+
+subroutine run_pt2(N_st,energy) 
+  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
+  !if(lambda_mrcc_pt2(0)  == 0) return
+  
+  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_mrcepa0_energy_pt2(energy(1)+pt2(1))
+
+end 
+
diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f
index 39b0f58e..f690c790 100644
--- a/src/Determinants/determinants.irp.f
+++ b/src/Determinants/determinants.irp.f
@@ -438,8 +438,12 @@ end
  do i=1,N_states
    psi_coef_min(i) = minval(psi_coef(:,i))
    psi_coef_max(i) = maxval(psi_coef(:,i))
-   abs_psi_coef_min(i) = dabs(psi_coef_min(i))
-   abs_psi_coef_max(i) = dabs(psi_coef_max(i))
+   abs_psi_coef_min(i) = minval( dabs(psi_coef(:,i)) )
+   abs_psi_coef_max(i) = maxval( dabs(psi_coef(:,i)) )
+   call write_double(6,psi_coef_max(i), 'Max coef')
+   call write_double(6,psi_coef_min(i), 'Min coef')
+   call write_double(6,abs_psi_coef_max(i), 'Max abs coef')
+   call write_double(6,abs_psi_coef_min(i), 'Min abs coef')
  enddo
 
 END_PROVIDER