diff --git a/.travis.yml b/.travis.yml
index f451f1d6..b5cf3053 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -24,7 +24,7 @@ python:
 
 script: 
   - ./configure  --production ./config/gfortran.cfg
-  - source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD
+  - source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD All_singles
   - source ./quantum_package.rc ; ninja
   - source ./quantum_package.rc ; cd ocaml ; make ; cd -
   - source ./quantum_package.rc ; cd tests ; bats bats/qp.bats
diff --git a/config/ifort.cfg b/config/ifort.cfg
index cc848cba..f8c36e4e 100644
--- a/config/ifort.cfg
+++ b/config/ifort.cfg
@@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32
 # 0 : Deactivate
 # 
 [OPTION]
-MODE    : OPT        ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
+MODE    : DEBUG      ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
 CACHE   : 1          ; Enable cache_compile.py
 OPENMP  : 1          ; Append OpenMP flags
 
diff --git a/plugins/Full_CI/full_ci.irp.f b/plugins/Full_CI/full_ci.irp.f
index 074ec7e1..56bf96d9 100644
--- a/plugins/Full_CI/full_ci.irp.f
+++ b/plugins/Full_CI/full_ci.irp.f
@@ -62,11 +62,27 @@ program full_ci
     endif
     print *,  'N_det          = ', N_det
     print *,  'N_states       = ', N_states
+    do  k = 1, N_states
+    print*,'State ',k
     print *,  'PT2            = ', pt2
     print *,  'E              = ', CI_energy
     print *,  'E(before)+PT2  = ', E_CI_before+pt2
+    enddo
     print *,  '-----'
     E_CI_before = CI_energy
+    if(N_states.gt.1)then
+     print*,'Variational Energy difference'
+     do i = 2, N_states
+      print*,'Delta E = ',CI_energy(i) - CI_energy(1)
+     enddo
+    endif
+    if(N_states.gt.1)then
+     print*,'Variational + perturbative Energy difference'
+     do i = 2, N_states
+      print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
+     enddo
+    endif
+    E_CI_before = CI_energy
     call ezfio_set_full_ci_energy(CI_energy)
     if (abort_all) then
       exit
diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg
index 4ab84b7a..b1c459ba 100644
--- a/src/Determinants/EZFIO.cfg
+++ b/src/Determinants/EZFIO.cfg
@@ -40,6 +40,12 @@ doc: Force the wave function to be an eigenfunction of S^2
 interface: ezfio,provider,ocaml
 default: False
 
+[diagonalize_s2]
+type: logical
+doc: Diagonalize the S^2 operator within the n_states_diag states required. Notice : the vectors are sorted by increasing S^2 values. 
+interface: ezfio,provider,ocaml
+default: True 
+
 [threshold_davidson]
 type: Threshold
 doc: Thresholds of Davidson's algorithm
diff --git a/src/Determinants/diagonalize_CI.irp.f b/src/Determinants/diagonalize_CI.irp.f
index 9af8d413..b533bed2 100644
--- a/src/Determinants/diagonalize_CI.irp.f
+++ b/src/Determinants/diagonalize_CI.irp.f
@@ -36,17 +36,36 @@ END_PROVIDER
  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
+  implicit none
+  double precision :: ovrlp,u_dot_v
+  integer :: i_good_state
+  integer, allocatable  :: index_good_state_array(:)
+  logical, allocatable  :: good_state_array(:)
+  double precision, allocatable :: s2_values_tmp(:)
+  integer :: i_other_state
+  double precision, allocatable  :: eigenvectors(:,:), eigenvalues(:)
+  integer :: i_state
+  double precision :: s2,e_0
+  integer                        :: i,j,k
+  double precision, allocatable :: s2_eigvalues(:)
+  double precision, allocatable :: e_array(:)
+  integer, allocatable :: iorder(:)
   
-  do j=1,N_states_diag
+  ! Guess values for the "N_states_diag" states of the CI_eigenvectors 
+  do j=1,min(N_states_diag,N_det)
     do i=1,N_det
       CI_eigenvectors(i,j) = psi_coef(i,j)
     enddo
   enddo
+
+  do j=N_det+1,N_states_diag
+    do i=1,N_det
+      CI_eigenvectors(i,j) = 0.d0
+    enddo
+  enddo
   
   if (diag_algorithm == "Davidson") then
     
@@ -59,36 +78,77 @@ END_PROVIDER
     
   else if (diag_algorithm == "Lapack") then
     
-    double precision, allocatable  :: eigenvectors(:,:), eigenvalues(:)
     allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
     allocate (eigenvalues(N_det))
     call lapack_diag(eigenvalues,eigenvectors,                       &
         H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
     CI_electronic_energy(:) = 0.d0
-    do i=1,N_det
-       CI_eigenvectors(i,1) = eigenvectors(i,1)
-    enddo
-    integer :: i_state
-    double precision :: s2
     if (s2_eig) then
       i_state = 0
+      allocate (s2_eigvalues(N_det))
+      allocate(index_good_state_array(N_det),good_state_array(N_det))
+      good_state_array = .False.
       do j=1,N_det
         call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2)
-        print*,'s2 = ',s2
+        s2_eigvalues(j) = s2
+        ! Select at least n_states states with S^2 values closed to "expected_s2"
         if(dabs(s2-expected_s2).le.0.3d0)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)
-        CI_eigenvectors_s2(i_state) = s2
+         i_state +=1
+         index_good_state_array(i_state) = j
+         good_state_array(j) = .True.
         endif
-        if (i_state.ge.N_states_diag) then
-          exit
+        if(i_state.eq.N_states) then
+         exit
         endif
       enddo
+      if(i_state .ne.0)then
+       ! Fill the first "i_state" states that have a correct S^2 value
+       do j = 1, i_state
+          do i=1,N_det
+           CI_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j))
+          enddo
+          CI_electronic_energy(j) = eigenvalues(index_good_state_array(j))
+          CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j))
+       enddo
+       i_other_state = 0
+       do j = 1, N_det 
+        if(good_state_array(j))cycle
+        i_other_state +=1
+        if(i_state+i_other_state.gt.n_states_diag)then
+         exit
+        endif
+        call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2)
+        do i=1,N_det
+         CI_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j)
+        enddo
+        CI_electronic_energy(i_state+i_other_state) = eigenvalues(j)
+        CI_eigenvectors_s2(i_state+i_other_state) = s2
+       enddo
+ 
+       deallocate(index_good_state_array,good_state_array)
+
+      else
+       print*,''
+       print*,'!!!!!!!!   WARNING  !!!!!!!!!'
+       print*,'  Within the ',N_det,'determinants selected'
+       print*,'  and the ',N_states_diag,'states requested'
+       print*,'  We did not find any state with S^2 values close to ',expected_s2
+       print*,'  We will then set the first N_states eigenvectors of the H matrix'
+       print*,'  as the CI_eigenvectors'
+       print*,'  You should consider more states and maybe ask for diagonalize_s2 to be .True. or just enlarge the CI space'
+       print*,''
+       do j=1,min(N_states_diag,N_det)
+         do i=1,N_det
+           CI_eigenvectors(i,j) = eigenvectors(i,j)
+         enddo
+        CI_electronic_energy(j) = eigenvalues(j)
+        CI_eigenvectors_s2(j) = s2_eigvalues(j)
+       enddo
+      endif
+     deallocate(s2_eigvalues)
     else
-      do j=1,N_states_diag
+      ! Select the "N_states_diag" states of lowest energy
+      do j=1,min(N_det,N_states_diag)
         call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
         do i=1,N_det
           CI_eigenvectors(i,j) = eigenvectors(i,j)
@@ -99,7 +159,100 @@ END_PROVIDER
     endif
     deallocate(eigenvectors,eigenvalues)
   endif
+
+  if(diagonalize_s2.and.n_states_diag > 1.and. n_det >= n_states_diag)then
+   ! Diagonalizing S^2 within the "n_states_diag" states found
+   allocate(s2_eigvalues(N_states_diag))
+   call diagonalize_s2_betweenstates(psi_det,CI_eigenvectors,n_det,size(psi_det,3),size(CI_eigenvectors,1),min(n_states_diag,n_det),s2_eigvalues)
+
+    do j = 1, N_states_diag
+     do i = 1, N_det
+      psi_coef(i,j) = CI_eigenvectors(i,j)
+     enddo
+    enddo
+
+   if(s2_eig)then
+
+   ! Browsing the "n_states_diag" states and getting the lowest in energy "n_states" ones that have the S^2 value 
+   ! closer to the "expected_s2" set as input
+
+    allocate(index_good_state_array(N_det),good_state_array(N_det))
+    good_state_array = .False.
+    i_state = 0
+    do j = 1, N_states_diag
+     if(dabs(s2_eigvalues(j)-expected_s2).le.0.3d0)then
+      good_state_array(j) = .True.
+      i_state +=1 
+      index_good_state_array(i_state) = j
+     endif
+    enddo
+    ! Sorting the i_state good states by energy 
+    allocate(e_array(i_state),iorder(i_state))
+    do j = 1, i_state
+     do i = 1, N_det 
+       CI_eigenvectors(i,j) = psi_coef(i,index_good_state_array(j))
+     enddo
+     CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j))
+     call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
+     CI_electronic_energy(j) = e_0
+     e_array(j) = e_0
+     iorder(j) = j
+    enddo
+    call dsort(e_array,iorder,i_state) 
+    do j = 1, i_state
+     CI_electronic_energy(j) = e_array(j)
+     CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(iorder(j)))
+     do i = 1, N_det
+       CI_eigenvectors(i,j) = psi_coef(i,index_good_state_array(iorder(j)))
+     enddo
+!    call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
+!    print*,'e    = ',CI_electronic_energy(j)
+!    print*,'<e>  = ',e_0
+!    call get_s2_u0(psi_det,CI_eigenvectors(1,j),N_det,size(CI_eigenvectors,1),s2)
+!    print*,'s^2  = ',CI_eigenvectors_s2(j)
+!    print*,'<s^2>= ',s2
+    enddo
+    deallocate(e_array,iorder)
+
+    ! Then setting the other states without any specific energy order
+    i_other_state = 0
+    do j = 1, N_states_diag
+     if(good_state_array(j))cycle
+     i_other_state +=1
+     do i = 1, N_det 
+       CI_eigenvectors(i,i_state + i_other_state) = psi_coef(i,j)
+     enddo
+     CI_eigenvectors_s2(i_state + i_other_state) = s2_eigvalues(j)
+     call u0_H_u_0(e_0,CI_eigenvectors(1,i_state + i_other_state),n_det,psi_det,N_int)
+     CI_electronic_energy(i_state + i_other_state) = e_0
+    enddo
+    deallocate(index_good_state_array,good_state_array)
+
+
+   else
+ 
+    ! Sorting the N_states_diag by energy, whatever the S^2 value is
+
+    allocate(e_array(n_states_diag),iorder(n_states_diag))
+    do j = 1, N_states_diag
+     call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
+     e_array(j) = e_0
+     iorder(j) = j
+    enddo
+    call dsort(e_array,iorder,n_states_diag) 
+    do j = 1, N_states_diag
+     CI_electronic_energy(j) = e_array(j)
+     do i = 1, N_det
+      CI_eigenvectors(i,j) = psi_coef(i,iorder(j))
+     enddo
+     CI_eigenvectors_s2(j) = s2_eigvalues(iorder(j))
+    enddo
+    deallocate(e_array,iorder)
+   endif
+   deallocate(s2_eigvalues)
+  endif
   
+ 
 END_PROVIDER
 
 subroutine diagonalize_CI
diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f
index 6da7b8ec..d84e4578 100644
--- a/src/Determinants/s2.irp.f
+++ b/src/Determinants/s2.irp.f
@@ -214,4 +214,175 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
   deallocate (shortcut, sort_idx, sorted, version)
 end
 
+subroutine get_uJ_s2_uI(psi_keys_tmp,psi_coefs_tmp,n,nmax_coefs,nmax_keys,s2,nstates)
+ implicit none
+ use bitmasks
+ integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax_keys)
+ integer, intent(in) :: n,nmax_coefs,nmax_keys,nstates
+ double precision, intent(in) :: psi_coefs_tmp(nmax_coefs,nstates)
+ double precision, intent(out) :: s2(nstates,nstates)
+ double precision :: s2_tmp,accu
+ integer :: i,j,l,jj,ll,kk
+ integer, allocatable           :: idx(:)
+ double precision, allocatable :: tmp(:,:)
+ BEGIN_DOC
+ ! returns the matrix elements of S^2 "s2(i,j)" between the "nstates" states 
+ ! psi_coefs_tmp(:,i) and psi_coefs_tmp(:,j)
+ END_DOC
+ s2 = 0.d0
+ do ll = 1, nstates
+  do jj = 1, nstates
+ accu = 0.d0
+ !$OMP PARALLEL DEFAULT(NONE)                                         &
+ !$OMP PRIVATE (i,j,kk,idx,tmp,s2_tmp) & 
+ !$OMP SHARED (ll,jj,psi_keys_tmp,psi_coefs_tmp,N_int,n,nstates)      &
+ !$OMP REDUCTION(+:accu)
+ allocate(idx(0:n))
+   !$OMP DO SCHEDULE(dynamic)
+   do i = 1, n
+    call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int)
+    accu += psi_coefs_tmp(i,ll) * s2_tmp * psi_coefs_tmp(i,jj)
+    call filter_connected(psi_keys_tmp,psi_keys_tmp(1,1,i),N_int,i-1,idx)
+    do kk=1,idx(0)
+     j = idx(kk)
+     call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int)
+     accu += psi_coefs_tmp(i,ll) * s2_tmp * psi_coefs_tmp(j,jj) + psi_coefs_tmp(i,jj) * s2_tmp * psi_coefs_tmp(j,ll)
+    enddo
+   enddo
+   !$OMP END DO NOWAIT
+ deallocate(idx)
+ !$OMP BARRIER
+ !$OMP END PARALLEL
+   s2(ll,jj) += accu
+  enddo
+ enddo
+ do i = 1, nstates
+  do j =i+1,nstates
+   accu = 0.5d0 * (s2(i,j) + s2(j,i))
+   s2(i,j) = accu
+   s2(j,i) = accu
+  enddo
+ enddo
+end
+
+subroutine diagonalize_s2_betweenstates(keys_tmp,psi_coefs_inout,n,nmax_keys,nmax_coefs,nstates,s2_eigvalues)
+ BEGIN_DOC
+! You enter with nstates vectors in psi_coefs_inout that may be coupled by S^2
+! The subroutine diagonalize the S^2 operator in the basis of these states. 
+! The vectors that you obtain in output are no more coupled by S^2, 
+! which does not necessary mean that they are eigenfunction of S^2. 
+! n,nmax,nstates = number of determinants, physical dimension of the arrays and number of states
+! keys_tmp = array of integer(bit_kind) that represents the determinants 
+! psi_coefs(i,j) = coeff of the ith determinant in the jth state
+! VECTORS ARE SUPPOSED TO BE ORTHONORMAL IN INPUT
+ END_DOC
+ implicit none
+ use bitmasks
+ integer, intent(in) :: n,nmax_keys,nmax_coefs,nstates
+ integer(bit_kind), intent(in) :: keys_tmp(N_int,2,nmax_keys)
+ double precision, intent(inout) :: psi_coefs_inout(nmax_coefs,nstates)
+
+!integer, intent(in) :: ndets_real,ndets_keys,ndets_coefs,nstates
+!integer(bit_kind), intent(in) :: keys_tmp(N_int,2,ndets_keys)
+!double precision, intent(inout) :: psi_coefs_inout(ndets_coefs,nstates)
+ double precision, intent(out)   :: s2_eigvalues(nstates)
+
+
+ double precision,allocatable :: s2(:,:),overlap(:,:)
+ double precision, allocatable :: eigvalues(:),eigvectors(:,:)
+ integer :: i,j,k
+ double precision, allocatable :: psi_coefs_tmp(:,:)
+ double precision :: accu,coef_contract
+ double precision :: u_dot_u,u_dot_v
+
+ print*,''
+ print*,'*********************************************************************'
+ print*,'Cleaning the various vectors by diagonalization of the S^2 matrix ...'
+ print*,''
+ print*,'nstates = ',nstates
+ allocate(s2(nstates,nstates),overlap(nstates,nstates))
+  do i = 1, nstates
+    overlap(i,i) = u_dot_u(psi_coefs_inout(1,i),n)
+    do j = i+1, nstates
+      overlap(i,j) = u_dot_v(psi_coefs_inout(1,j),psi_coefs_inout(1,i),n)
+      overlap(j,i) = overlap(i,j)
+    enddo
+  enddo
+ print*,'Overlap matrix in the basis of the states considered'
+ do i = 1, nstates
+  write(*,'(10(F16.10,X))')overlap(i,:)
+ enddo
+ call ortho_lowdin(overlap,size(overlap,1),nstates,psi_coefs_inout,size(psi_coefs_inout,1),n)
+ print*,'passed ortho'
+
+  do i = 1, nstates
+    overlap(i,i) = u_dot_u(psi_coefs_inout(1,i),n)
+    do j = i+1, nstates
+      overlap(i,j) = u_dot_v(psi_coefs_inout(1,j),psi_coefs_inout(1,i),n)
+      overlap(j,i) = overlap(i,j)
+    enddo
+  enddo
+ print*,'Overlap matrix in the basis of the Lowdin orthonormalized states '
+ do i = 1, nstates
+  write(*,'(10(F16.10,X))')overlap(i,:)
+ enddo
+
+ call get_uJ_s2_uI(keys_tmp,psi_coefs_inout,n_det,size(psi_coefs_inout,1),size(keys_tmp,3),s2,nstates)
+ print*,'S^2 matrix in the basis of the states considered'
+ double precision :: accu_precision_diag,accu_precision_of_diag
+ accu_precision_diag = 0.d0
+ accu_precision_of_diag = 0.d0
+ do i = 1, nstates
+  do j = i+1, nstates
+   if(  ( dabs(s2(i,i) - s2(j,j)) .le.1.d-10 ) .and. (dabs(s2(i,j) + dabs(s2(i,j)))) .le.1.d-10) then
+    s2(i,j) = 0.d0
+    s2(j,i) = 0.d0
+   endif
+  enddo
+ enddo
+ do i = 1, nstates
+  write(*,'(10(F10.6,X))')s2(i,:)
+ enddo
+
+ print*,'Diagonalizing the S^2 matrix'
+
+ allocate(eigvalues(nstates),eigvectors(nstates,nstates))
+ call lapack_diagd(eigvalues,eigvectors,s2,nstates,nstates)
+ print*,'Eigenvalues of s^2'
+ do i = 1, nstates
+  print*,'s2 = ',eigvalues(i)
+  s2_eigvalues(i) = eigvalues(i)
+ enddo
+
+ print*,'Building the eigenvectors of the S^2 matrix'
+ allocate(psi_coefs_tmp(nmax_coefs,nstates))
+ psi_coefs_tmp = 0.d0
+ do j = 1, nstates
+  do k = 1, nstates
+   coef_contract =  eigvectors(k,j)    !  <phi_k|Psi_j>
+   do i = 1, n_det
+    psi_coefs_tmp(i,j) += psi_coefs_inout(i,k) * coef_contract
+   enddo
+  enddo
+ enddo
+ do j = 1, nstates
+  accu = 0.d0
+   do i = 1, n_det
+    accu += psi_coefs_tmp(i,j) * psi_coefs_tmp(i,j)
+   enddo
+   print*,'Norm of vector = ',accu
+   accu = 1.d0/dsqrt(accu)
+   do i = 1, n_det
+    psi_coefs_inout(i,j) = psi_coefs_tmp(i,j) * accu
+   enddo
+ enddo
+!call get_uJ_s2_uI(keys_tmp,psi_coefs_inout,n_det,size(psi_coefs_inout,1),size(keys_tmp,3),s2,nstates)
+!print*,'S^2 matrix in the basis of the NEW states considered'
+!do i = 1, nstates
+! write(*,'(10(F16.10,X))')s2(i,:)
+!enddo
+
+ deallocate(s2,eigvalues,eigvectors,psi_coefs_tmp,overlap)
+
+end
 
diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f
index 202d310d..e983ec34 100644
--- a/src/Determinants/slater_rules.irp.f
+++ b/src/Determinants/slater_rules.irp.f
@@ -1507,6 +1507,33 @@ subroutine get_occ_from_key(key,occ,Nint)
   
 end
 
+subroutine u0_H_u_0(e_0,u_0,n,keys_tmp,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Computes e_0 = <u_0|H|u_0>/<u_0|u_0>
+  !
+  ! n : number of determinants
+  !
+  END_DOC
+  integer, intent(in)            :: n,Nint
+  double precision, intent(out)  :: e_0
+  double precision, intent(in)   :: u_0(n)
+  integer(bit_kind),intent(in)   :: keys_tmp(Nint,2,n)
+  
+  double precision               :: H_jj(n)
+  double precision               :: v_0(n)
+  double precision               :: u_dot_u,u_dot_v,diag_H_mat_elem
+  integer :: i,j
+  do i = 1, n
+   H_jj(i) = diag_H_mat_elem(keys_tmp(1,1,i),Nint)
+  enddo
+  
+  call H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
+  e_0 = u_dot_v(v_0,u_0,n)/u_dot_u(u_0,n)
+end
+
+
 subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
   use bitmasks
   implicit none