diff --git a/data/ezfio_defaults b/data/ezfio_defaults
index 1de27eff..be7e7605 100644
--- a/data/ezfio_defaults
+++ b/data/ezfio_defaults
@@ -17,8 +17,9 @@ cis_dressed
 
 determinants
   n_states                   1
+  n_states_diag              determinants_n_states
   n_det_max_jacobi           5000
-  threshold_generators       0.999
+  threshold_generators       0.995
   threshold_selectors        0.999
   read_wf                    False
 
diff --git a/src/Dets/SC2.irp.f b/src/Dets/SC2.irp.f
index aa2f49d1..921cf896 100644
--- a/src/Dets/SC2.irp.f
+++ b/src/Dets/SC2.irp.f
@@ -175,7 +175,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence)
      enddo
      call lapack_diag(eigenvalues,eigenvectors,                       &
          H_matrix_tmp,size(H_matrix_all_dets,1),sze)
-     do j=1,min(N_states,sze)
+     do j=1,min(N_states_diag,sze)
        do i=1,sze
          u_in(i,j) = eigenvectors(i,j)
        enddo
diff --git a/src/Dets/davidson.irp.f b/src/Dets/davidson.irp.f
index a79e9123..1efb70ff 100644
--- a/src/Dets/davidson.irp.f
+++ b/src/Dets/davidson.irp.f
@@ -12,7 +12,7 @@ BEGIN_PROVIDER [ integer, davidson_sze_max ]
   ! Max number of Davidson sizes
   END_DOC
   ASSERT (davidson_sze_max <= davidson_iter_max)
-  davidson_sze_max = 8*N_states
+  davidson_sze_max = 8*N_states_diag
 END_PROVIDER
 
 subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit)
diff --git a/src/Dets/determinants.ezfio_config b/src/Dets/determinants.ezfio_config
index ff97ca36..4641db0f 100644
--- a/src/Dets/determinants.ezfio_config
+++ b/src/Dets/determinants.ezfio_config
@@ -4,6 +4,7 @@ determinants
   mo_label      character*(64)
   n_det         integer
   n_states      integer
+  n_states_diag      integer
   psi_coef      double precision (determinants_n_det,determinants_n_states)
   psi_det       integer*8 (determinants_n_int*determinants_bit_kind/8,2,determinants_n_det)
   n_det_max_jacobi integer
diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f
index f7ec2dfc..1b52fafd 100644
--- a/src/Dets/determinants.irp.f
+++ b/src/Dets/determinants.irp.f
@@ -150,7 +150,7 @@ subroutine read_dets(det,Nint,Ndet)
 end
 
 
-BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
+BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states_diag) ]
   implicit none
   BEGIN_DOC
   ! The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file
@@ -161,6 +161,11 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
   logical                        :: exists
   double precision, allocatable  :: psi_coef_read(:,:)
   character*(64)                 :: label
+
+  psi_coef = 0.d0
+  do i=1,N_states_diag
+    psi_coef(i,i) = 1.d0
+  enddo
   
   if (read_wf) then
     call ezfio_has_determinants_psi_coef(exists)
@@ -183,22 +188,8 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
       enddo
       deallocate(psi_coef_read)
       
-    else
-      
-      psi_coef = 0.d0
-      do i=1,N_states
-        psi_coef(i,i) = 1.d0
-      enddo
-
     endif
     
-  else
-    
-    psi_coef = 0.d0
-    do i=1,N_states
-      psi_coef(i,i) = 1.d0
-    enddo
-    
   endif
     
   
diff --git a/src/Dets/diagonalize_CI.irp.f b/src/Dets/diagonalize_CI.irp.f
index 1329b581..f3a352f3 100644
--- a/src/Dets/diagonalize_CI.irp.f
+++ b/src/Dets/diagonalize_CI.irp.f
@@ -9,13 +9,13 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ]
     diag_algorithm = "Lapack"
   endif
 
-  if (N_det < N_states) then
+  if (N_det < N_states_diag) then
     diag_algorithm = "Lapack"
   endif
   
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ]
+BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ]
   implicit none
   BEGIN_DOC
   ! N_states lowest eigenvalues of the CI matrix
@@ -24,23 +24,25 @@ BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ]
   integer                        :: j
   character*(8)                  :: st
   call write_time(output_Dets)
-  do j=1,N_states
+  do j=1,N_states_diag
     CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion
     write(st,'(I4)') j
     call write_double(output_Dets,CI_energy(j),'Energy of state '//trim(st))
+    call write_double(output_Dets,CI_eigenvectors_s2(j),'S^2 of state '//trim(st))
   enddo
 
 END_PROVIDER
 
- BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states) ]
-&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states) ]
+ BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states_diag) ]
+&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states_diag) ]
+&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2, (N_states_diag) ]
   implicit none
   BEGIN_DOC
   ! Eigenvectors/values of the CI matrix
   END_DOC
   integer                        :: i,j
   
-  do j=1,N_states
+  do j=1,N_states_diag
     do i=1,N_det
       CI_eigenvectors(i,j) = psi_coef(i,j)
     enddo
@@ -49,7 +51,7 @@ END_PROVIDER
   if (diag_algorithm == "Davidson") then
     
     call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy, &
-        size(CI_eigenvectors,1),N_det,N_states,N_int,output_Dets)
+        size(CI_eigenvectors,1),N_det,N_states_diag,N_int,output_Dets)
     
   else if (diag_algorithm == "Lapack") then
     
@@ -63,17 +65,18 @@ END_PROVIDER
     double precision :: s2
     j=0
     i_state = 0
-!   do while(i_state.lt.N_states)
+    do while(i_state.lt.min(N_states_diag,N_det))
       j+=1
-!     call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
-!     if(dabs(s2-expected_s2).le.0.1d0)then
+      call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
+      if(dabs(s2-expected_s2).le.0.1d0)then
        i_state += 1
        do i=1,N_det
          CI_eigenvectors(i,i_state) = eigenvectors(i,j)
        enddo
        CI_electronic_energy(i_state) = eigenvalues(j)
-!     endif
-!   enddo
+       CI_eigenvectors_s2(i_state) = s2
+      endif
+    enddo
     deallocate(eigenvectors,eigenvalues)
   endif
   
@@ -86,10 +89,10 @@ subroutine diagonalize_CI
 !  eigenstates of the CI matrix
   END_DOC
   integer :: i,j
-  do j=1,N_states
+  do j=1,N_states_diag
     do i=1,N_det
       psi_coef(i,j) = CI_eigenvectors(i,j)
     enddo
   enddo
-  SOFT_TOUCH psi_coef CI_electronic_energy CI_energy CI_eigenvectors
+  SOFT_TOUCH psi_coef CI_electronic_energy CI_energy CI_eigenvectors CI_eigenvectors_s2
 end
diff --git a/src/Dets/diagonalize_CI_SC2.irp.f b/src/Dets/diagonalize_CI_SC2.irp.f
index 412dc2ae..02252826 100644
--- a/src/Dets/diagonalize_CI_SC2.irp.f
+++ b/src/Dets/diagonalize_CI_SC2.irp.f
@@ -1,13 +1,13 @@
-BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states) ]
+BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states_diag) ]
   implicit none
   BEGIN_DOC
-  ! N_states lowest eigenvalues of the CI matrix
+  ! N_states_diag lowest eigenvalues of the CI matrix
   END_DOC
   
   integer                        :: j
   character*(8)                  :: st
   call write_time(output_Dets)
-  do j=1,N_states
+  do j=1,N_states_diag
     CI_SC2_energy(j) = CI_SC2_electronic_energy(j) + nuclear_repulsion
     write(st,'(I4)') j
     call write_double(output_Dets,CI_SC2_energy(j),'Energy of state '//trim(st))
@@ -23,15 +23,15 @@ END_PROVIDER
  threshold_convergence_SC2 = 1.d-10
 
  END_PROVIDER
- BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states) ]
-&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states) ]
+ BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ]
+&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states_diag) ]
   implicit none
   BEGIN_DOC
   ! Eigenvectors/values of the CI matrix
   END_DOC
   integer                        :: i,j
   
-  do j=1,N_states
+  do j=1,N_states_diag
     do i=1,N_det
       CI_SC2_eigenvectors(i,j) = psi_coef(i,j)
 !     CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j)
@@ -41,17 +41,17 @@ END_PROVIDER
   
   double precision :: convergence
   call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, &
-        size(CI_SC2_eigenvectors,1),N_det,N_states,N_int,threshold_convergence_SC2)
+        size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
 END_PROVIDER
 
 subroutine diagonalize_CI_SC2
   implicit none
   BEGIN_DOC
-!  Replace the coefficients of the CI states by the coefficients of the 
+!  Replace the coefficients of the CI states_diag by the coefficients of the 
 !  eigenstates of the CI matrix
   END_DOC
   integer :: i,j
-  do j=1,N_states
+  do j=1,N_states_diag
     do i=1,N_det
       psi_coef(i,j) = CI_SC2_eigenvectors(i,j)
     enddo
diff --git a/src/Dets/options.irp.f b/src/Dets/options.irp.f
index df2ba1b5..c918e2c2 100644
--- a/src/Dets/options.irp.f
+++ b/src/Dets/options.irp.f
@@ -9,6 +9,16 @@ T.set_ezfio_name( "N_states" )
 T.set_output    ( "output_dets" )
 print T
 
+
+T.set_type      ( "integer" )
+T.set_name      ( "N_states_diag" )
+T.set_doc       ( "Number of states to consider for the diagonalization " )
+T.set_ezfio_dir ( "determinants" )
+T.set_ezfio_name( "N_states_diag" )
+T.set_output    ( "output_dets" )
+print T
+
+
 T.set_name      ( "N_det_max_jacobi" )
 T.set_doc       ( "Maximum number of determinants diagonalized by Jacobi" )
 T.set_ezfio_name( "N_det_max_jacobi" )
diff --git a/src/Dets/s2.irp.f b/src/Dets/s2.irp.f
index 6993675c..18c3dfaf 100644
--- a/src/Dets/s2.irp.f
+++ b/src/Dets/s2.irp.f
@@ -52,7 +52,7 @@ BEGIN_PROVIDER [ double precision, expected_s2]
      call ezfio_get_determinants_expected_s2(expected_s2)
    else
      expected_s2 = elec_alpha_num - elec_beta_num + 0.5d0 * ((elec_alpha_num - elec_beta_num)**2*0.5d0 - (elec_alpha_num-elec_beta_num))
-     call ezfio_set_determinants_expected_s2(expected_s2)
+!    call ezfio_set_determinants_expected_s2(expected_s2)
    endif
 
 END_PROVIDER