diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py
index 8b947259..53a44dae 100755
--- a/scripts/generate_h_apply.py
+++ b/scripts/generate_h_apply.py
@@ -141,7 +141,7 @@ class H_apply(object):
       """
       self.data["size_max"] = "256" 
       self.data["initialization"] = """
-      PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors E_corr_per_selectors
+      PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
       """
       self.data["keys_work"] = """
       call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
@@ -206,7 +206,6 @@ class H_apply(object):
       call copy_h_apply_buffer_to_wf
       selection_criterion_min = selection_criterion_min*0.1d0
       selection_criterion = selection_criterion_min
-      !call remove_small_contributions
       """
       self.data["keys_work"] = """
       e_2_pert_buffer = 0.d0
diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f
index d36efb9b..5c6ce54a 100644
--- a/src/Dets/H_apply_template.f
+++ b/src/Dets/H_apply_template.f
@@ -341,7 +341,7 @@ subroutine $subroutine($params_main)
   
   $decls_main
   
-  PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map N_det_reference psi_generators
+  PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map N_det_selectors psi_generators
   integer                        :: i_generator, k
   double precision               :: wall_0, wall_1, wall_2
   
diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f
index 31b82136..aa7a7a31 100644
--- a/src/Dets/connected_to_ref.irp.f
+++ b/src/Dets/connected_to_ref.irp.f
@@ -1,10 +1,77 @@
-integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
+logical function is_in_wavefunction(key,keys,Nint,N_past_in,Ndet)
+  implicit none
+
+  integer, intent(in)            :: Nint, N_past_in, Ndet
+  integer(bit_kind), intent(in)  :: keys(Nint,2,Ndet)
+  integer(bit_kind), intent(in)  :: key(Nint,2)
+
+  integer                        :: i, ibegin, iend, istep, l
+  integer*8                      :: det_ref, det_search
+  integer*8, external            :: det_search_key
+  
+  is_in_wavefunction = .False.
+  ibegin = 1
+  iend   = N_det
+  ASSERT (N_past > 0)
+  ASSERT (N_det >= N_past)
+  
+  det_ref = det_search_key(key,Nint)
+  
+  istep = ishft(iend-ibegin+1,-1)
+  i=ibegin+istep
+  do while (istep > 1)
+    i = ibegin + istep
+    det_search = det_search_key(psi_det_sorted_bit(1,1,i),Nint)
+!   print *,  istep, det_ref, det_search
+    if ( det_search > det_ref ) then
+      iend = i
+    else if ( det_search == det_ref ) then
+      exit
+    else
+      ibegin = i
+    endif
+    istep = ishft(iend-ibegin+1,-1)
+    i = ibegin + istep
+  end do
+! pause
+
+  do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref)
+    i = i-1
+    if (i == 0) then
+      exit
+    endif
+  enddo
+  i += 1
+  do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref)
+    if ( (key(1,1) /= psi_det_sorted_bit(1,1,i)).or.                               &
+          (key(1,2) /= psi_det_sorted_bit(1,2,i)) ) then
+      continue
+    else
+      is_in_wavefunction = .True.
+      !DEC$ LOOP COUNT MIN(3)
+      do l=2,Nint
+        if ( (key(l,1) /= psi_det_sorted_bit(l,1,i)).or.                           &
+              (key(l,2) /= psi_det_sorted_bit(l,2,i)) ) then
+          is_in_wavefunction = .False.
+          exit
+        endif
+      enddo
+      if (is_in_wavefunction) then
+        return
+      endif
+    endif
+    i += 1
+    
+  enddo
+  
+end
+
+integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
   use bitmasks
   implicit none
   integer, intent(in)            :: Nint, N_past_in, Ndet
   integer(bit_kind), intent(in)  :: keys(Nint,2,Ndet)
   integer(bit_kind), intent(in)  :: key(Nint,2)
-  double precision, intent(in)   :: thresh
   
   integer                        :: N_past
   integer                        :: i, l
@@ -44,15 +111,21 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
       return
     endif
     
-    do i=N_past,Ndet
-      if ( (key(1,1) /= keys(1,1,i)).or.                             &
-            (key(1,2) /= keys(1,2,i)) ) then
-        cycle
-      endif
-      connected_to_ref = -i
-      return
-    enddo
+    logical, external :: is_in_wavefunction
+    if (is_in_wavefunction(key,keys,Nint,N_past_in,Ndet)) then
+      connected_to_ref = -1
+    endif
     return
+
+  ! do i=N_past,Ndet
+  !   if ( (key(1,1) /= keys(1,1,i)).or.                             &
+  !         (key(1,2) /= keys(1,2,i)) ) then
+  !     cycle
+  !   endif
+  !   connected_to_ref = -i
+  !   return
+  ! enddo
+  ! return
     
   else if (Nint==2) then
     
@@ -80,7 +153,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
     endif
     
     !DIR$ LOOP COUNT (1000)
-    do i=N_past+1,Ndet
+    do i=N_past,Ndet
       if ( (key(1,1) /= keys(1,1,i)).or.                             &
             (key(1,2) /= keys(1,2,i)).or.                            &
             (key(2,1) /= keys(2,1,i)).or.                            &
@@ -119,7 +192,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
       return
     endif
     
-    do i=N_past+1,Ndet
+    do i=N_past,Ndet
       if ( (key(1,1) /= keys(1,1,i)).or.                             &
             (key(1,2) /= keys(1,2,i)).or.                            &
             (key(2,1) /= keys(2,1,i)).or.                            &
@@ -161,7 +234,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
       return
     endif
     
-    do i=N_past+1,Ndet
+    do i=N_past,Ndet
       if ( (key(1,1) /= keys(1,1,i)).or.                             &
             (key(1,2) /= keys(1,2,i)) ) then
         cycle
diff --git a/src/Dets/determinants.ezfio_config b/src/Dets/determinants.ezfio_config
index f4b8b908..a6e2f4c2 100644
--- a/src/Dets/determinants.ezfio_config
+++ b/src/Dets/determinants.ezfio_config
@@ -7,4 +7,6 @@ determinants
   psi_det       integer*8 (determinants_N_int*determinants_bit_kind/8,2,determinants_n_det)
   H_apply_threshold  double precision
   n_det_max_jacobi integer
+  threshold_generators  double precision
+  threshold_selectors   double precision
 
diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f
index 45f0c3f0..570741cc 100644
--- a/src/Dets/determinants.irp.f
+++ b/src/Dets/determinants.irp.f
@@ -90,15 +90,6 @@ END_PROVIDER
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ integer, N_det_reference ]
- implicit none
- BEGIN_DOC
- ! Number of determinants in the reference wave function
- END_DOC
-   N_det_reference = N_det
- ASSERT (N_det_reference > 0)
-END_PROVIDER
-
 BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (N_det) ]
  implicit none
  BEGIN_DOC
@@ -123,7 +114,7 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted, (N_det) ]
  implicit none
  BEGIN_DOC
- ! Wave function sorted by determinants (state-averaged)
+ ! Wave function sorted by determinants contribution to the norm (state-averaged)
  END_DOC
  integer :: i,j,k
  integer, allocatable ::  iorder(:)
@@ -149,3 +140,54 @@ END_PROVIDER
 
 END_PROVIDER
 
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_bit, (N_int,2,N_det) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_sorted_bit, (N_det,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Determinants on which we apply <i|H|psi> for perturbation.
+ !o They are sorted by determinants interpreted as integers. Useful
+ ! to accelerate the search of a determinant
+ END_DOC
+ integer :: i,j,k
+ integer, allocatable ::  iorder(:)
+ integer*8, allocatable :: bit_tmp(:)
+ integer*8, external :: det_search_key
+
+ allocate ( iorder(N_det), bit_tmp(N_det) )
+
+ do i=1,N_det
+   iorder(i) = i
+   !$DIR FORCEINLINE
+   bit_tmp(i) = det_search_key(psi_det(1,1,i),N_int)
+ enddo
+ call i8sort(bit_tmp,iorder,N_det)
+ !DIR$ IVDEP
+ do i=1,N_det
+  do j=1,N_int
+    psi_det_sorted_bit(j,1,i) = psi_det(j,1,iorder(i))
+    psi_det_sorted_bit(j,2,i) = psi_det(j,2,iorder(i))
+  enddo
+  do k=1,N_states
+    psi_coef_sorted_bit(i,k) = psi_coef(iorder(i),k)
+  enddo
+ enddo
+
+ deallocate(iorder)
+
+END_PROVIDER
+
+
+integer*8 function det_search_key(det,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Return an integer*8 corresponding to a determinant index for searching
+  END_DOC
+  integer, intent(in) :: Nint
+  integer(bit_kind), intent(in) :: det(Nint,2)
+  integer :: i
+  det_search_key = iand(det(1,1),det(1,2))
+  do i=2,Nint
+    det_search_key = ieor(det_search_key,iand(det(i,1),det(i,2)))
+  enddo
+end
diff --git a/src/Full_CI/full_ci.irp.f b/src/Full_CI/full_ci.irp.f
index 72eec40e..3651ecea 100644
--- a/src/Full_CI/full_ci.irp.f
+++ b/src/Full_CI/full_ci.irp.f
@@ -24,5 +24,14 @@ program cisd
       exit
     endif
   enddo
+  print *,  'Final step'
+  call remove_small_contributions
+  call diagonalize_CI
+  print *,  'N_det    = ', N_det
+  print *,  'N_states = ', N_states
+  print *,  'PT2      = ', pt2
+  print *,  'E        = ', CI_energy
+  print *,  'E+PT2    = ', CI_energy+pt2
+  print *,  '-----'
   deallocate(pt2,norm_pert)
 end
diff --git a/src/Generators_full/generators.irp.f b/src/Generators_full/generators.irp.f
index 00930ae6..40013131 100644
--- a/src/Generators_full/generators.irp.f
+++ b/src/Generators_full/generators.irp.f
@@ -1,14 +1,25 @@
 use bitmasks
 
-BEGIN_PROVIDER [ double precision, generators_threshold ]
+BEGIN_PROVIDER [ double precision, threshold_generators ]
  implicit none
  BEGIN_DOC
  ! Percentage of the norm of the state-averaged wave function to
  ! consider for the generators
  END_DOC
- generators_threshold = 0.97d0
+ logical                        :: exists
+ PROVIDE ezfio_filename
+ call ezfio_has_determinants_threshold_generators(exists)
+ if (exists) then
+   call ezfio_get_determinants_threshold_generators(threshold_generators)
+ else
+   threshold_generators = 0.99d0
+   call ezfio_set_determinants_threshold_generators(threshold_generators)
+ endif
+ ASSERT (N_det > 0)
+ call write_double(output_Dets,threshold_generators,'Threshold on generators')
 END_PROVIDER
 
+
 BEGIN_PROVIDER [ integer, N_det_generators ]
  implicit none
  BEGIN_DOC
@@ -22,7 +33,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
  N_det_generators = N_det
  do i=1,N_det
    norm = norm + psi_average_norm_contrib_sorted(i)
-   if (norm > generators_threshold) then
+   if (norm > threshold_generators) then
      N_det_generators = i-1
      exit
    endif
diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES
index 5cd66cf8..0ef6daa9 100644
--- a/src/NEEDED_MODULES
+++ b/src/NEEDED_MODULES
@@ -1 +1,2 @@
-AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation SingleRefMethod CISD_selected Selectors_full MP2 Generators_full Full_CI CISD_SC2_selected CISD_SC2_selected 
+AOs BiInts Bitmask CISD CISD_SC2_selected CISD_selected DensityMatrix Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Perturbation SC2 Selectors_full SingleRefMethod Utils
+
diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f
index d96b557a..b73053ef 100644
--- a/src/Perturbation/perturbation_template.f
+++ b/src/Perturbation/perturbation_template.f
@@ -15,7 +15,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
   double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
   double precision               :: c_pert(N_st), e_2_pert(N_st),  H_pert_diag(N_st)
   integer                        :: i,k, c_ref
-  integer                        :: connected_to_ref
+  integer, external              :: connected_to_ref
   
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
@@ -24,7 +24,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
   ASSERT (N_st > 0)
   do i = 1,buffer_size
 
-    c_ref = connected_to_ref(buffer(1,1,i),psi_generators,Nint,i_generator,N_det,h_apply_threshold)
+    c_ref = connected_to_ref(buffer(1,1,i),psi_generators,Nint,i_generator,N_det)
 
     if (c_ref /= 0) then
       cycle
diff --git a/src/Selectors_full/selectors.irp.f b/src/Selectors_full/selectors.irp.f
index 8f963f0d..e430f591 100644
--- a/src/Selectors_full/selectors.irp.f
+++ b/src/Selectors_full/selectors.irp.f
@@ -1,31 +1,68 @@
 use bitmasks
 
+BEGIN_PROVIDER [ double precision, threshold_selectors ]
+ implicit none
+ BEGIN_DOC
+ ! Percentage of the norm of the state-averaged wave function to
+ ! consider for the selectors
+ END_DOC
+ logical                        :: exists
+ PROVIDE ezfio_filename
+ call ezfio_has_determinants_threshold_selectors(exists)
+ if (exists) then
+   call ezfio_get_determinants_threshold_selectors(threshold_selectors)
+ else
+   threshold_selectors = 0.99d0
+   call ezfio_set_determinants_threshold_selectors(threshold_selectors)
+ endif
+ ASSERT (N_det > 0)
+ call write_double(output_Dets,threshold_selectors,'Threshold on selectors')
+END_PROVIDER
+
 BEGIN_PROVIDER [ integer, psi_selectors_size ]
  implicit none
  psi_selectors_size = psi_det_size
 END_PROVIDER
+
 BEGIN_PROVIDER [ integer, N_det_selectors]
  implicit none
-  N_det_selectors = N_det
+ BEGIN_DOC
+ ! For Single reference wave functions, the number of selectors is 1 : the
+ ! Hartree-Fock determinant
+ END_DOC
+ integer :: i
+ double precision :: norm
+ call write_time(output_dets)
+ norm = 0.d0
+ N_det_selectors = N_det
+ do i=1,N_det
+   norm = norm + psi_average_norm_contrib_sorted(i)
+   if (norm > threshold_selectors) then
+     N_det_selectors = i-1
+     exit
+   endif
+ enddo
+ N_det_selectors = max(N_det_selectors,1)
+ call write_int(output_dets,N_det_selectors,'Number of selectors')
 END_PROVIDER
 
  BEGIN_PROVIDER [ integer(bit_kind), psi_selectors, (N_int,2,psi_selectors_size) ]
 &BEGIN_PROVIDER [ double precision, psi_selectors_coef, (psi_selectors_size,N_states) ]
   implicit none
   BEGIN_DOC
-  ! On what we apply <i|H|psi> for perturbation. If selection, it may be 0.9 of the norm.
+  ! Determinants on which we apply <i|H|psi> for perturbation.
   END_DOC
   integer                        :: i,k
 
   do i=1,N_det_selectors
     do k=1,N_int
-      psi_selectors(k,1,i) = psi_det(k,1,i)
-      psi_selectors(k,2,i) = psi_det(k,2,i)
+      psi_selectors(k,1,i) = psi_det_sorted(k,1,i)
+      psi_selectors(k,2,i) = psi_det_sorted(k,2,i)
     enddo
   enddo
   do k=1,N_states
     do i=1,N_det_selectors
-      psi_selectors_coef(i,k) = psi_coef(i,k)
+      psi_selectors_coef(i,k) = psi_coef_sorted(i,k)
     enddo
   enddo
 END_PROVIDER