diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f
index a5e85777..aa748628 100644
--- a/src/davidson/diagonalization_hs2_dressed.irp.f
+++ b/src/davidson/diagonalization_hs2_dressed.irp.f
@@ -450,7 +450,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       if (s2_eig) then
         h_p = s_
         do k=1,shift2
-          h_p(k,k) = h_p(k,k) + S_z2_Sz - expected_s2
+          h_p(k,k) = h_p(k,k) - expected_s2
         enddo
         if (only_expected_s2) then
           alpha = 0.1d0
@@ -496,7 +496,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
           0.d0, s_, size(s_,1))
 
       do k=1,shift2
-        s2(k) = s_(k,k) + S_z2_Sz
+        s2(k) = s_(k,k)
       enddo
 
       if (only_expected_s2) then
diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f
index 996011cd..33e9478a 100644
--- a/src/davidson/diagonalize_ci.irp.f
+++ b/src/davidson/diagonalize_ci.irp.f
@@ -107,7 +107,7 @@ END_PROVIDER
        H_prime(1:N_det,1:N_det) = H_matrix_all_dets(1:N_det,1:N_det) +  &
          alpha * S2_matrix_all_dets(1:N_det,1:N_det)
        do j=1,N_det
-         H_prime(j,j) = H_prime(j,j) + alpha*(S_z2_Sz - expected_s2)
+         H_prime(j,j) = H_prime(j,j) - alpha*expected_s2
        enddo
        call lapack_diag(eigenvalues,eigenvectors,H_prime,size(H_prime,1),N_det)
        CI_electronic_energy(:) = 0.d0
diff --git a/src/determinants/s2.irp.f b/src/determinants/s2.irp.f
index 391d0073..d73b2dbf 100644
--- a/src/determinants/s2.irp.f
+++ b/src/determinants/s2.irp.f
@@ -8,24 +8,35 @@ double precision function diag_S_mat_elem(key_i,Nint)
   BEGIN_DOC
 ! Returns <i|S^2|i>
   END_DOC
-  integer                        :: nup, i
-  integer(bit_kind)              :: xorvec(N_int_max)
+  integer                        :: nup, ntot, i
+  integer(bit_kind)              :: xorvec(N_int_max), upvec(N_int_max)
 
   do i=1,Nint
     xorvec(i) = xor(key_i(i,1),key_i(i,2))
   enddo
 
   do i=1,Nint
-    xorvec(i) = iand(xorvec(i),key_i(i,1))
+    upvec(i) = iand(xorvec(i),key_i(i,1))
   enddo
 
+  ! nup is number of alpha unpaired
+  ! ntot is total number of unpaired
   nup = 0
+  ntot = 0
   do i=1,Nint
     if (xorvec(i) /= 0_bit_kind) then
-      nup += popcnt(xorvec(i))
+      ntot += popcnt(xorvec(i))
+      if (upvec(i) /= 0_bit_kind) then
+        nup += popcnt(upvec(i))
+      endif
     endif
   enddo
-  diag_S_mat_elem = dble(nup)
+  
+  double precision :: sz
+  sz = nup - 0.5d0*ntot
+  
+  !<S^2> = <S+ S-> + Sz(Sz-1)
+  diag_S_mat_elem = nup + sz*(sz-1)
 
 end
 
@@ -125,7 +136,7 @@ subroutine u_0_S2_u_0(e_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
 
   call S2_u_0_nstates(v_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
   do i=1,N_st
-    e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n) + S_z2_Sz
+    e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n)
   enddo
 end